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Preface 



Wind tunnels are the primary research tools used in aerodynamic research. They are 
used to study the effects of air moving past solid objects. Although great advances in 
computational methods have been made in recent years, wind tunnel tests remain es- 
sential for obtaining the full range of data required to guide detailed design decisions 
for various practical engineering problems. 

This book collects original and innovative research studies on recent applications in 
wind tunnel tests, exhibiting various investigation directions and providing a bird's 
eye view on this broad subject area. It is composed of seven chapters that have been 
grouped in two major parts. The first part of the book (chapters 1-4) deals with wind 
tunnel technologies and devices. The second part (chapters 5-7) deals with the latest 
applications of wind tunnel testing. 

The following is a brief description of the subjects that are covered in each chapter: 

Chapter 1 reviews some examples of environmental wind tunnels. 
Chapter 2 describes a 6-DOF system for the measurements of forces and torques in 
wind tunnels. 

Chapter 3 proposes a 6-DOF wire-driven parallel manipulator with redundant actua- 
tions for wind tunnels. 

Chapter 4 introduces the plasma wind tunnel test on a large thermal protection system 
demonstrator. 

Chapter 5 describes the flow visualization and the proper orthogonal decomposition 
of aeroelastic phenomena. 

Chapter 6 introduces the wind tunnel testing of pneumatic artificial muscles. 
Chapter 7 provides the flow-induced vibrations and scattering of roof tiles by wind 
tunnel testing. 

The text is addressed not only to researchers but also to professional engineers, engi- 
neering lecturers, and students seeking to gain better understanding of the current 
status of wind tunnels. 

Through its seven chapters, the reader will have an access to a wide range of works 
related to wind tunnel testing. 
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Environmental Wind Tunnels 

Jonathan Merrison 

Aarhus University, 
Denmark 



1. Introduction 



Wind tunnels have been used extensively in industry and research applications over the 
past 50 years. They vary greatly in scale and geometry, with some large enough to house 
and test small aircraft (see for example NASA, ATP facilities) and others are miniaturized 
flow generators used in the calibration of small sensors. However they invariably utilize the 
same basic technology and design elements. Similarly environmental simulators are also 
used widely in research, for example in climate and planetary studies. Here again they 
superficially vary greatly in size and configuration, but basically consist of a hermetic 
chamber with some form of temperature control [Jensen et al. 2008]. There is therefore a 
broad array of standard and often commercial technologies and construction techniques 
which have been successfully applied within the fields of wind tunnel and environmental 
simulator design. Some of these technologies and techniques will be outlined in this chapter 
to aid researchers or technology developers in their efforts to design or use environmental 
wind tunnels and also serve as an informative guide to those new to these fields of 
investigation. 

The fusion of an environmental simulator and a wind tunnel is a natural evolution of 
laboratory based technology to fulfill the need to reproduce specific physical conditions 
found in nature. Although facilities of this kind are only now being fully developed, they 
have the potential to expand into a new research field that could substantially contribute to 
our understanding of climate and mediate growth in advanced sensor technologies. In this 
chapter many of the challenges in designing and constructing environmental wind tunnels 
will be introduced and possible solutions presented, with some emphasis placed on extreme 
terrestrial and Martian planetary conditions. In addition some of the many and varied 
scientific and industrial applications will be discussed. Generally environmental wind 
tunnels are already in current use as a method of testing and calibrating meteorology 
sensors of various kinds especially wind flow sensors (anemometers). Application of wind 
tunnels in civil engineering and town planning is becoming common place. Here through 
wind tunnel simulation and modeling the flow of air around buildings and through built-up 
areas may be useful to avoid the generation of high wind shear and hazardous vortices at 
periods of high wind or storms. Such simulations can also aid in the design and placement 
of wind generation systems such as wind turbines. 

The formalized scaling laws developed by Reynolds (Reynolds equations) allows 
measurements, for example in smaller scale laboratory wind tunnels, which generate the 
same (or extremely similar) flow to that generated in the natural setting [Monin and Yaglom 
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1973, Hall 1988, Mollinger and Nieuwstadt 1996, Fay and Sonwalkar 1991]. This scaling law 
involves the relationship between wind speed, spatial scale and viscosity such that adjusting 
and combining these parameters can allow realistic laboratory simulation for example on 
the cm-m scale of flow dynamics on the 10s to 100s of meters. It can also, for example, allow 
comparison of effects in one fluid (e.g. air) to be translated into those seen in another fluid 
such as water. This technique has been successfully applied in the design of all forms of 
transport, such as aircraft, ships and cars. 

Wind tunnel studies have and are contributing powerfully in attempts to understand and 
describe the action of wind in arid areas. Following the pioneering work of Bagnold, including 
the use of laboratory (and field) wind tunnels, the study of Aeolian (wind driven) sand 
transport has evolved into a scientific research field [Bagnold 1941]. It is now clear that 
Aeolian transport has a great impact on local environments and on the global climate through 
the production of aerosols, the erosion of surface material and the serious environmental 
problem of desertification. Aeolian transport of sand/ dust under planetary conditions other 
than Earths is also of great importance to understanding these extreme environments and can 
help achieve a deeper understanding of our own environment. For example Aeolian processes 
are seen on Mars, Venus and Saturn's moon Titan, but are probably found on any planetary 
body with a significant atmosphere. Sand features such as dunes are common on these planets 
and in the case of Mars dust entrainment is seen to be the most powerful climatic factor. 
Interesting differences in the Aeolian features seen in these extra-terrestrial environments is 
the spatial scale compared to those on Earth. The study of extra terrestrial Aeolian phenomena 
can only effectively be studied in the laboratory using an environmental wind tunnel 
simulator. Even with such simulators, only some aspects of Aeolian transport on other planets 
can be successfully reproduced, such as the surface shear stress, wind speed, fluid density, 
temperature, humidity and (more ambitiously) surface microstructure, adhesive properties. 
Other physical aspects are extremely problematic, for example gravity and specifics of the 
surface composition (mineralogy). 

An obvious application for an environmental wind tunnel is the study of the upper 
atmosphere (the troposphere and stratosphere), specifically low temperatures, low 
pressures and the presence of aerosols of various types. Clearly this is of relevance to the 
aircraft industry, especially (high altitude) jet aircraft. The recent (2010) disturbance in 
Atlantic flights due to the generation of dust aerosols by the Icelandic volcano 
(Eyjafjallajokull) is a good example, where a deeper understanding of these aerosols in the 
upper atmosphere could possibly have avoided a large degree of disruption. The 
development of new aerosol sensor technologies also appears to be necessary. In fact wind 
tunnels can both help to unravel the complex dynamics of aerosol behavior and to 
understand their formation processes through the generation of fine suspended mineral 
particulates (dust). It should be stressed here that the study of aerosols is far from being 
limited to a global climatic factor. Aerosols present a real hazard to environmental and 
human safety both in the home and in local environments. Conversely aerosols are also used 
widely in medicine and the pharmaceutical and cosmetics industries. Specifically nano- 
micro meter scale particulates suspended in the air can penetrate the deep lung as well as be 
suspended for long periods of time (months) in the atmosphere and transported great 
distances (globally). Smoke, clouds, dust, are just some of the many forms of aerosol that 
affect our environment and can be studied in environmental wind tunnels to better 
understand their (apparently complex) behavior as well as develop new technology in order 
to quantify and control them. 
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Fig. 1. Left upper and lower; satellite photographs of Mars and Earth, North Africa 
respectively, showing dust storms and clouds. (Courtesy NASA/ JPL-Cal tech). Right upper 
fog in Mars, Valles Marineris taken by the High Resolution Stereo Camera (HRSC) on board 
ESA's Mars Express spacecraft. Right lower; acid haze seen amongst the thick clouds of 
Venus, photographed by the ESA's Venus Express spacecraft. 

In the future the study of aerosols will probably be the single most important application of 
environmental wind tunnels and it is hoped that the work presented here will contribute 
towards these types of study. 



2. Environmental wind tunnel mechanical design 

There are two basic types of wind tunnel design which may be referred to as Open Circuit 
or Closed Circuit (or closed cycle). In a terrestrial (ambient pressure) open circuit wind 
tunnel design fresh air is drawn (or blown) into the entrance and expelled at the exit, 
whereas in a closed circuit wind tunnel the expelled air is fed again into the inlet such that 
the same air is re-circulated. Either of these two wind tunnel types can be housed in an 
environmental (or planetary simulation) chamber giving rise to two distinct types of 
environmental wind tunnel design. These different wind tunnel types (shown schematically 
in figures 2a-2d) have distinct characteristics, their advantages and disadvantages will be 
discussed. 

The implementation of thermal and flow control within these differing system designs will 
vary. In the case of the open circuit design flow and thermal control systems should be 
implemented upwind and focus primarily on manipulating the gas which is inlet. In the 
case of a re-circulating design, since the system is a closed cycle, flow correction and thermal 
control can in principle be implemented in any (or all) sections of the circuit. In practice the 
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implementation of thermal control will depend on the thermal control system chosen and 
general technical restraints of the wind tunnel design. Similarly flow control will depend on 
the desired flow characteristics and the practical limitations on resources. 
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Fig. 2. Different types of Wind Tunnel geometry combining open/ closed circuit designs and 
ambient or enclosed (environmental control). 

In traditional ambient pressure wind tunnel facilities the choice of construction materials is 
largely unrestricted. Materials are therefore chosen dependent on mechanical properties 
(strength, weight, etc.) and possibly also cost and availability, wood for example is used in 
many wind tunnels. For environmental wind tunnels the choice of materials is generally far 
more restrictive since, to maintain low pressure or gas purity, materials with low out- 
gassing properties should be chosen and for temperature control the thermal properties and 
mechanical properties at low temperatures must be considered. The choice of materials 
subsequently affects the mechanical design of the wind tunnel structure. 
External access to the wind tunnel (especially the test section) is also of great importance in 
most cases, both during operation and installation or maintenance. Here access includes, for 
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Fig. 3. Schematics of Left; Aarhus University Wind Tunnel I (AWTSI) design. Center; AWTS 
II design. Right; open circuit ambient showing upwind flow control (see flow control) 

example mechanical, electrical and optical (visual) systems. Specifically mechanical access 
could involve being able to orientate a sample or sensor and therefore require rotation or 
translation mechanisms. Electrical access may be in the form of cabling for power and data 
transfer. Optical access could be cameras, lighting, spectrometers or other optical sensors. 
Ideally these forms of access should be as spatially close to the active section of the wind 
tunnel as possible and preferably large in cross section. For ideal flow (i.e. minimizing 
boundary effects) a wind tunnel should be cylindrical in cross section, however for the 
housing and access of samples/ sensors, as well as many other practical applications of wind 
tunnels, it is desirable to use a rectangular cross section. This does not constitute a problem 
for most ambient-pressure applications; however for an enclosed (pressurized) wind tunnel 
this does present a technical challenge. The two environmental wind tunnel systems at 
Aarhus University apply two radically different geometrical solutions to this problem, with 
the AWTS-I system housing the cylindrical wind tunnel within its own (cylindrical) return 
flow, giving a rather attractive flow transport and uniform cross-section, though poor access 
to the test section and a non-optimal (circular) cross-section [Merrison et al. 2008]. The 
AWTS-II design conversely has an attractive, almost rectangular wind tunnel cross-section 
and good access to the test section, however the return flow is divided into two, above and 
below the test section, giving extremely non-ideal flow and constriction of the flow in the 
return section, which resulted in the need for extensive flow correction. At low pressure 
(below lOOmbar) the highest wind speed achieved by AWTSI is around 15-20m/s, whereas 
AWTSII has achieved 20-25m/s, with similar degrees of (free flow) turbulence for both wind 
tunnels i.e. 5-20% increasing with wind speed. 




Fig. 4. Photographs of the (10m long) AWTS-II facility showing the mobile environmental 
chamber sections, the central test section can be removed laterally. 

In contrast to the Aarhus environmental wind tunnels the NASA Ames MARSWIT (Mars 
Surface Wind Tunnel, California USA) is an open-circuit, low pressure wind tunnel 
powered by a high pressure nozzle ejector system, the total length is 13m with a main test 
section of 1.2m by 0.9 m and is housed in a 4000 m^ low-pressure chamber which can 
operate at pressures down to -3.8 mbar and wind speeds of 20m/ s - 180m/ s (at low 
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pressure) [White 1981, Greeley and Iversen 1985]. This system cannot be cooled and has 
been used for boundary layer studies. 

For low pressure wind tunnel systems the structure of the vacuum chamber is one of the 
primary design features. This will typically require the use of a thick (bulky) steel shell and 
frame which, for mechanical strength, will optimally be cylindrical/ spherical in form. This 
is similarly true for high pressure vessels. For open circuit environmental wind tunnels the 
limitations on the pressure vessel will limit the size and geometry of the test section. 
However for a re-circulating environmental wind tunnel the pressure vessel will even more 
strongly restrict design of the wind tunnel since, assuming the largest free flow cross section 
is desired then there must still be sufficient space for the return flow to be housed. It is 
desirable for this return flow cross section to be comparable to the test section cross section 
to avoid high turbulence and turbulent losses. The AWTS-II facility is one of the largest 
environmental wind tunnels with a cross section of around 2mxlm and a chamber volume 
of around 40m3, it is significantly larger than the almost Im^ volume and cross section of 
0.4mx0.4m of the AWTS-I. 




Fig. 5. A Light Emitting Diode based light source for solar simulation, illumination or crude 
spectroscopy in an environmental wind tunnel. Upper Left shows a photograph of a single 
array section including one of each of the seven different (wavelength) LEDs, Upper Right 
shows the irradiance measured within the wind tunnel (with all LEDs activated) showing 
the single wavelength components, the Lower photographs are taken inside the wind tunnel 
test section as different colored LED arrays are activated (red, green blue), the LED array 
strips are mounted in the upper two edges of this section. 



In both industrial and scientific applications a common requirement is a light source which 
simulates the solar irradiance over a broad wavelength range. A problem with many light 
sources, for example halogen lamps and discharge lamps, is the generation of heat, both 
conductive and as thermal radiation (infra-red) which can make environmental temperature 
control difficult. Employing a light source outside the environmental chamber alleviates this 
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problem, however it then restricts the illumination of samples considerably. Compromise 
here will generally be necessary. An attractive option is the use of light emitting diodes 
(LEDs) which are efficient and monochromatic, being available as intense sources though 
generating relatively little heat. LEDs are low voltage making them technically easy to 
implement in most cases. With the use of arrays of variously colored LED the correct light 
irradiance can be achieved within broad optical wavelengths and even into the near infra- 
red (more than lOOOnm) and the Ultra Violet, with the latest UV LEDs below 250nm. 

3. Flow control 

Wind tunnels vary in their requirements for flow uniformity, while some are designed for 
low turbulence (close to laminar) flow, others apply techniques in order to reproduce 
particular boundary layer conditions (often referred to as boundary layer wind tunnels) 
whereas for some a specific free-flow degree of turbulence is required. In these differing 
cases it is probably fair to say that they are attempting to reproduce differing turbulence 
regimes present in nature and that it is therefore difficult to generalize about these wind 
flow designs. However it is worth discussing differing flow control techniques which are 
commonly employed and how specifically they can be applied. 

Flow guides are smooth plates of differing geometry which are installed in order to steer the 
wind flow to obtain a desired wind pattern. For example they may be; curved in order to 
guide the flow around bends, they may be planar in order to straighten the flow or they may 
be used to partition the flow into sections in order to prevent unwanted lateral flow/ eddies. 
In open circuit wind tunnels flow control should (obviously) be installed upwind. However, 
in a re-circulating wind tunnel they should generally be installed at the source of the 
unwanted flow pattern, which could be upwind or down wind. In the case of the European 
Mars Environmental Wind Tunnel (see figure 6) flow guides have been used to great effect 
at the entrance to the wind generating fans system and prevented extremely destructive 
back-flow caused by the rotation of the fan blades. Meshes are used to reduce turbulence in 
the wind flow and to obtain a more homogeneous flow profile, especially on scales larger 
than the mesh size which is typically of the order of 1mm. This is done at the expense of 
wind speed. Meshes are often utilized as a set of two separated by some mm-cm. In this case 
a pressure gradient is generated across the meshes, this helps to disrupt turbulence and non- 
uniformities in the flow. It should be noted that both flow guides and meshes while 
improving flow properties, will typically increase friction and therefore reduce the (net) 
wind flow for a particular wind generation power. The use of upwind roughness blocks and 
turbulence spires manipulate the vertical wind flow profile (at the test section) in order to 
emulate an infinite upwind 'fetch' i.e. to reproduce the surface boundary layer flow which 
would be produced if the wind tunnel were infinite in length. Clearly this is of great 
importance when studying boundary layer effects such as the entrainment and transport of 
sand or the flow patterns around a surface feature [Irwin 1981, Shao and Raupach 1992]. 
Expansion and compression stages can be used in wind tunnel design to increase wind 
speed, improve flow linearity and reduce turbulence. Here compression of the wind tunnel 
will increase the downwind flow speed and reduce the relative transverse turbulence. 
Clearly this is done at the cost of wind tunnel cross-sectional size and is not always possible 
to implement especially within a re-circulating wind tunnel. Often in open circuit wind 
tunnels and invariably in re-circulating systems wind generation is provided by a fan or 
fans. Fan design is in many cases non trivial, involving modeling and calculation regarding 
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the specific choice of fan blade size, number, form, angle and also motor power, torque and 
rotation rate. Such modeling and calculation can be aided by computational fluid dynamic 
calculations. Here one begins with the required parameters of wind speed (and ambient 
pressure), based on the wind tunnel design. The flow calculations will then predict a certain 
degree of frictional loss as a function of wind speed. The fan system can then be modeled as 
a system to generate a pressure gradient necessary to balance this frictional loss and 
maintain the desired flow rate. Given the flow rate and the required pressure differential a 
particular fan design can be chosen i.e. these are the required input parameters for the 
choice of fan design. In the case of environmental wind tunnels the choice of fan material 
must also be considered, for example to be compatible with out-gassing limits and low/ high 
temperature. 
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Fig. 6. Photographs into the flow generation section at the AWTS-II facility. Left shows the 
1.8m diameter fans installed. Center shows with the upper and lower flow separators and 
Right the system of vertical and horizontal flow guides compartmentalizing the flow, 
preventing rotation and excessive turbulence. 

Since almost all forms of high power motor are incompatible with the demands of (low) 
pressure and temperature within an environmental chamber, the drive mechanism for a fan 
system must be mounted externally. This presents a problem for the transfer of torque to the 
fan since passing a rapidly rotating axel through a pressure seal system is also incompatible 
with avoiding pressure leaks and maintaining low temperatures. A possible solution which 
has been employed in the various facilities at Aarhus University is the use of a magnetic 
coupling [Merrison et al. 2008]. Such couplings are commercially available and transfer 
torque from the drive axel (external) to the fan axel (internal) through a complex of magnetic 
fields generated by permanent magnets. This avoids physical contact of the two axels and 
allows this coupling to be completely hermetic (vacuum tight). A drawback with this system 
is the limited degree of torque which can be transferred by such couplings before they begin 
to slip which may limit the rotation rate (wind speed) within the wind tunnel. It does 
however have the benefit of protecting the drive-fan system from damage as slippage of this 
coupling is not hazardous. 

A type of open circuit environmental chamber has been employed for Mars simulation 
conditions at Oxford University. Here gas is injected from an array of (relatively small) 
inlets into a flow volume which is continually being evacuated by a pump. In this case an 
extremely low turbulence flow can be achieved along with high flow speeds as well as 
cooling. A drawback can be that the flow rate is dependent upon the chamber pressure such 
that control of low flow speed involves inlet and pump rate control. Such a system can be 
well suited to anemometer calibration and high wind speed tests [Wilson et al. 2008]. 
Discussion here has focused on low wind speeds (subsonic flows). There are however, forms 
of wind tunnel which generate and utilize supersonic and even hyper sonic flows for 
various studies. Specific applications are in the design and testing of supersonic aircraft or 



Environmental Wind Tunnels 



11 



re-entry devices. It should be noted that such wind tunnels utilize specialized techniques 
and the flow in such high velocity regimes differs from that at wind speeds significantly 
below that of sound [Barlow 1999]. Generally environmental wind tunnels will involve 
compromising the 'ideal' wind flow characteristics due to geometric constraints imposed by 
the environmental chamber or environmental control systems, for example reduced cross 
section, increased turbulence, reduced maximum wind speed or the use of cumbersome 
flow control systems. 




Fig. 7. Computational Fluid Dynamic calculations of an object within a wind tunnel 
showing: Upper; the finite element structure. Center; the calculated wind speed flow from 
red (high) to blue (low) and Lower; suspended (aerosol) particulates added to the flow and 
their trajectories traced. 



4. Computational fluid dynamics 

This chapter has focused upon experimental/ laboratory studies using environmental wind 
tunnels, however discussion should be made of the use of computational fluid dynamic 
modeling in this regard as in some cases this may be an alternative to laboratory simulation. 
However in most cases these two techniques are complementary. When constructing a fluid 
dynamic model in order to perform computational flow analysis it is necessary to make 
simplifications and assumptions which in most cases must be verified experimentally in 
order for confidence to be placed on the results [Peric et al. 1999]. A specific example is the 
calculation of flow around an irregular shaped object. In this case it is necessary to construct 
a finite element representation of this geometry before inputting wind flow boundary 
conditions. Although the resolution of this finite element array can be increased in order to 
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ascertain convergence, this will also be limited by computing power. Here comparison with 
experiment can be of great benefit in identifying sources of high sensitivity in the flow such 
that resolution be enhanced in this volume (see figure 7). A combination of targeted 
laboratory measurements and computational analysis can be ideal in simulating complex 
and difficult flow problems [Kinch et al. 2005]. Typically CFD is employed in the design 
phase of wind tunnels, though often the flow is complex and multi-dimensional such that 
empirical measurement and the implementation of correction elements is necessary to arrive 
at the most satisfactory flow characteristics. 

5. Environmental sensing technology 

A crucial aspect to any application of wind tunnels and/ or environmental simulators is the 
use of accurate and reliable sensor systems for control and reproducibility of the simulated 
conditions. Some sensor systems are readily and commercially available at a well evolved 
level, for example for temperature and pressure. Other sensor systems can be complex, 
expensive and require adaptation, examples are wind sensors (anemometers) and gas 
composition. For some sensor systems there is a clear demand for new technology, yet this 
technology awaits development, examples are shear stress sensing and aerosol analysis. 
In temperature sensing ther mo-resistors are widely available (for example 100 Ohm 
platinum resistors i.e. PtlOO), these are typically inexpensive and are accurate (typically 
around 1°C) over a wide range. The same could also be said of thermocouples (e.g. K-type), 
though these generally have a limited low temperature range. Thermocouples can also be 
difficult to integrate into an environmental chamber due to the need to maintain the contact 
potential i.e. maintain the exotic metal cables. Pressure sensor systems are available either 
for high pressure use, low pressures or specific to terrestrial conditions, i.e. limited to 
around 1 bar. Low pressure sensors are typically (generically) referred to as vacuum gauges. 
A type of vacuum gauge which is ideal for moderate low pressures (down to say O.lmbar) 
and which is accurate even in differing gas compositions is the capacitance vacuum sensor, 
it is therefore well suited to study of Earth's upper atmosphere or Mars. Such capacitor 
based techniques are also useful for determining pressure differentials which can be 
important in wind tunnel design or wind sensing (see Pi tot tube). Although absolute 
humidity (water vapor pressure) sensors are typically complex and expensive, relative 
humidity sensors are often extremely compact and operate over wide temperature and 
pressure ranges. Specifically thin polymer film type sensors are commercially available and 
are easily implemented into an environmental system (e.g. Honeywell HIH series). 
In environmental systems where the atmospheric composition may be controlled it is 
important to be able to monitor it. There are few available options in this case and typically a 
sensor system called a Rest Gas Analyzer is used. These are often a type of quadrapole (radio 
frequency) mass spectrometer. They operate by ionizing the gas at low pressure (i.e. leaked 
through a valve) and extracting the ion fragments individually to determine their mass to 
charge ratio. It may then be possible to re-construct the original molecular structure of the gas, 
it is however difficult if the atmosphere contains several species where some fragments are 
ambiguous and it is often difficult to precisely determine abundances without careful 
control/ calibration of the system and some expertise. Although these systems are relatively 
expensive and cumbersome to install, there is at present a lack of viable alternatives. 
Finally in any application where an array of sensory systems is used, it is desirable to 
implement a data-logging system which records the various sensor outputs during 
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measurement cycles. For environmental wind tunnel systems it is also natural then to 
integrate this data logging capability into a computer system which also interfaces (and 
records) some of the control parameters of the facility such as wind generation (driving fan 
rotation rate), vacuum/ pressure control system (pumps, valves etc,), cooling/ heating 
systems or lighting subsystems. Although this constitutes an added level of complexity it 
allows for a higher level of reproducibility, sensor correlation and possibly safety. 

6. Flow sensing technology 

Clearly of primary importance with regard to wind tunnels is the accurate sensing of wind 
flow (Anemometry). There is a wide variety of available anemometer techniques, some 
dating back over 500 years, others are still being developed. These wind sensing systems 
vary in accuracy, complexity, price, size, and so on. In the following paragraphs some of the 
most common wind sensing technologies will be presented and briefly discussed, 
specifically with respect to their application in wind tunnels. 




Fig. 8. Photographs of Laser Anemometers (acting also as suspended dust sensors). Left 
prototype time of flight instrument. Center the sensor during aerosol testing in an 
environmental wind tunnel. Right A commercial Laser Doppler Anemometer operating 
through an environmental wind tunnel access window, note the beams illuminating the 
suspended dust in the flow. 



6.1 Laser anemometers 

These are probably the most advanced and desirable type of wind sensor which have been 
applied in wind tunnels, specifically the Laser Doppler Anemometer (LDA) is used 
extensively. Several more recent variations on this instrument can measure in multiple 
dimensions, image and determine suspended grain size. This technique has the benefit of 
being non contact, such that it is independent of the environmental conditions within the 
flow (pressure, temperature, composition, etc.), it is also accurate and does not normally 
require external calibration. In fact LDA based systems are widely used in wind tunnel 
applications for the calibration of other types of wind sensor. The principle behind the 
technique is the scattering and detection of light by suspended aerosol particles, by 
measuring the frequency shift due to the velocity induced Doppler effect. More specifically 
two (or more) beams are use to produce an interference pattern, measurement of the shift in 
this pattern allows single velocity components of the grains to be determined. The system 
does have the disadvantage of requiring the presence of suspended particulates within the 
flow, which are introduced as smoke in many systems. However, for systems studying 
aerosols this is a major advantage since the suspended grain concentration can be quantified 
using this technique. Typically LDA systems are expensive and bulky, though can use 
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optical fibers and therefore achieve a relatively compact sensing head. Miniature (even 
micro-scale) laser based wind sensors are being developed, though have yet to advance 
from prototyping. One such system is based on a time of flight principle in which a light 
pattern is generated within the sensing volume. Single suspended aerosol particulates 
traversing this light pattern will scatter light with a modulated signal from which its 
velocity can be established, specifically in the case of the prototype shown in figure 8 a three 
line light pattern is used and the scattered light signal will consist of three pulses the time 
separation is then directly related to the velocity [Merrison et al. 2004, Merrison et al. 2006]. 
This type of technology has the potential to become miniaturized (on the sub-cm scale) and 
have low power consumption as well as being robust. Although limited in precision 
compared to LDA systems, it may be applied in systems too small or inaccessible for larger 
sensors and provide an affordable (and portable/ battery driven) aerosol sensor. The current 
advancements in solid state laser and other optoelectronic technology give sensors of this 
kind a promising future. 

6.2 Mechanical (cup anemometers or wind socks) 

Mechanical anemometers are by the far the oldest, simplest, most common and varied form 
of wind sensor. Most widely used are cup anemometers and forms of wind sock or wind 
vane. A cup anemometer consists typically of conical cups mounted on a axel such that 
wind drag causes rotational motion which can be sensed by a tachometer in order to relate 
the rotation rate to the wind speed. Wind vanes and socks are typically more primitive and 
consist of a structure (tube/ sock or plate) which is deflected by the wind such that the 
deflection angle is a measure of the wind speed and the direction may often be seen in the 
direction of the deflection. Such mechanical wind sensors are rarely used in wind tunnel 
applications due to their poor accuracy/ precision and often limited dynamic range. They 
are however an invariable component of weather/ climatic stations on Earth and have even 
been adapted for the extreme environment of Mars and Venus. Such systems can potentially 
be extremely compact, light weight, sensitive and robust given careful design and testing 
[Gunnlaugsson et al. 2008]. 

6.3 Hot wire or hot film 

These sensors have been used extensively in wind tunnel experiments over several decades. 
They are typically accurate and sensitive in terrestrial conditions, they can also be multi 
dimensional and have reasonably fast response times. Compared to mechanical wind 
sensing techniques they therefore provide improvement in precision. The measurement 
technique relies on (electrically) heating a thin wire or foil which is then cooled by the flow 
of air. The cooling rate is therefore related to the wind speed. There are many variations on 
the this concept including specialized geometries, multiple heated elements (to determine 
wind direction), pulsed operation and heater-sensor feedback circuitry. Challenges to this 
technique are thermal (conductive) losses and temperature dependences in addition to the 
sensitivity to atmospheric properties. Also the heated sensors are often physically fragile 
and poorly suited to harsh environments. However it has been demonstrated that careful 
design, testing and importantly calibration can allow these sensors to be used even in low 
pressure, thermally unstable environments such as Mars. The first successful wind sensor 
system developed by NASA was such a hot film anemometer. 



Environmental Wind Tunnels 15 

6.4 Pitot tubes 

Pitot tubes are a simple and widely applied wind velocity sensor. This type of sensor is used 
in the aerospace industry (airplanes) as well as wind tunnels. The principle is measuring the 
overpressure generated in a wind facing tube compared to a non wind facing aperture. This 
pressure differential is a function of the wind speed relative to the tube. It is therefore well 
suited to situations where the direction of the wind flow is known. Despite their wide use, 
the Pitot tube is typically limited in range (due to its strong dependence upon wind speed) 
and requires careful calibration, since it is dependent upon atmospheric conditions 
(pressure, temperature, etc.). 

6.5 Sonic anemometers 

Sonic anemometers are a relatively modern and commercially available sensor for determining 
wind flow, they utilize the transmission of high frequency sound (ultrasonic) in order to 
measure wind flow by determining the acoustic propagation speed. Sonic anemometers can 
simultaneously measure wind velocity in all three dimensions and at high sampling rate. 
These sensors are precise and being three dimensional are capable of quantifying vertical as 
well as lateral flow rates. This makes them the instrument of choice for the study of boundary 
layer transport. They are currently used widely in climatic/ atmospheric studies, though not 
usually in wind tunnel applications. Unfortunately sonic anemometers are sensitive to the 
physical properties of the atmosphere (composition, pressure, temperature, humidity etc.). 
This makes them poorly suited to many environmental applications. Research groups have 
attempted to adapt sonic anemometers to extreme environments such as that on Mars, though 
have been hindered by the low pressure. 

6.6 Shear stress 

The quantification of surface shear stress within a wind tunnel is crucially important when 
trying to evaluate the threshold or transport rates of granular material or more generally 
mass transport rates or heat transfer. Currently a large body of semi-empirical work allows 
the measurement of surface wind velocity to be related to the surface shear stress (friction 
velocity). More crudely measurement of the wind velocity, turbulence and surface 
roughness can be used to obtain estimates of shear stress [White 1991]. However 
experimentally these are often difficult and indirect approaches to the determination of 
surface shear stress. Ideally the application of nano-micro scale force/ pressure sensors could 
now allow the direct measurement of wind shear stress [Xu et al 2003], however these are 
not commercially available and have not advanced from research prototypes. 

7. Thermal control 

Most of the discussion here will concern cooling within environmental wind tunnels rather 
than heating, though in many respects the problems and solutions are essentially the same. 
In industry environmental wind tunnels typically refer to wind tunnels within which the 
temperature can be controlled, with heating and cooling over the range typically expected 
on earth i.e. around -60°C to +50°C, though with no control of pressure. Such wind tunnels 
are used extensively in the automobile and aerospace industries and are often on a scale 
(many square meters cross section) such that full size vehicles can be housed. In this case 
commercial refrigeration (freezer) technology can be employed. Cooling systems vary 
depending on the temperature range and power requirements, typically for temperatures 
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above -80°C closed cycle refrigeration (using refrigerants such as haloalkanes, ammonia, 
alcohol, etc.) can be used. Below this and down to around -190°C it is common to use liquid 
nitrogen flow through systems of some kind, since this is readily available and relatively 
cheap [Jensen et al. 2008]. Cryogenic pumps (e.g. closed cycle helium refrigeration systems) 
are used for very low temperatures (a few K) or smaller systems which operate for long 
periods of time and require stability. A convenient cooling or rather heat exchange system is 
the Peltier element. This is a low voltage DC driven thermo-electric device for generating a 
temperature gradient across thin plate type bimetallic networks. Thermal power exchange 
of the order of lOOW and temperature gradients of the order of 70°C (from room 
temperature) are typical. These are easily installed and relatively cheap and often used for 
thermal control of small detectors, solid state lasers, etc. A limitation is that they operate 
poorly at low temperature with a decrease in operating voltage (power) and temperature 
gradient, for example at around -130°C the maximum temperature gradient possible falls to 
around 10°C. In the case of heating systems, although electrical heaters are commonly 
available and a well established technology, their use in an environmental simulator 
imposes certain restrictions on the applicable materials and geometry, for example that the 
material should not outgas, should be compatible with cryogenic temperatures, be fully 
electrically encapsulated and available in various geometries depending on the wind tunnel 
structure. A commercially available solution which has many of these qualities are custom 
manufactured encapsulated heater mats (X-Mat®) they consist of copper conductors 
encased in Capton foil (3 mm thick) with silicone encased cables. 

In addition to an effective cooling/ heating system there is the need for efficient transport of 
thermal power to the sample or test section from the cooling system and possibly also to the 
ambient wind tunnel gas. Similarly there will be the desire to prevent thermal transport to 
the sample/ test section from external sources or even heating due to power generated by 
electrical systems. Thermal transport is therefore of importance and an understanding of it 
necessary. Three thermal transport processes are relevant in wind tunnel and environmental 
chamber design. Conductive thermal transport is generally of most importance and the 
method by which most thermal control can be achieved. Specifically, metals with high 
thermal conductivity are used, such as copper or aluminum. When dealing with fluids (such 
as gas) thermal transport can also occur by convective transport, here thermal power is 
transported by the physical flow of mass. This must be considered even when dealing with 
low pressure gases (a few mbar) since thermal conductivity is not strongly dependent upon 
pressure/ mass density and is therefore still significant at low pressure. An added 
complication is the possibility for volatiles, such as water, to condense/ freeze, evaporate 
and transfer heat through this mechanism. Radiative thermal transport i.e. the transport of 
heat through light emission/ absorption (typically at infrared wavelengths) is significant 
even in the absence of fluid/ solid contact. It can be of comparable importance in the thermal 
balance of samples or gas within the wind tunnel. Manipulation of this effect can be 
achieved by alteration of the surface properties of materials within the environmental 
chamber for example with the use of coatings (paints). 

For extreme temperatures (especially cryogenic) environmental chambers are necessarily 
going to have to employ some form of thermal insulation, preferably with good efficiency 
such that there is as little thermal contact (loss) between the test section and the ambient 
environment as possible. There are many types of ambient pressure insulation material 
available, varying in insulation efficiency, price, mass, volume etc. For low pressure systems 
the use of double walled (vacuum) insulation techniques (such as that used in Dewar flasks) 
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is effective and widely applied, however this technique essentially requires two nested 
vacuum systems and is therefore generally expensive and complex to construct. Efficient 
vacuum compatible thermal (cryogenic) insulation is available. One example is multi-layer 
thin film super-insulation developed at CERN-CryoLab and commercially available 
(JEHIER). Despite its high price relative to ambient pressure insulation, it is simple to apply 
multilayer super insulation within an environmental chamber design, it is also reasonably 
efficient and affordable compared to a multi-walled vacuum insulation solution. 
In a closed circuit systems the frictional loss of power within wind tunnel can cause 
significant heating at elevated wind speeds (of the order lOOW/m^). This can become 
problematic for thermal control systems in such situations. Specifically heating of the gas 
will ensue and heat deposition on surfaces in contact with the gas. This must be considered 
when designing thermal control systems if stable temperatures are to be achieved. 
As mentioned previously with respect to sensors, in addition to a thermal sensor system and 
a cooling/ heating system an automated (intelligent) control network needs to support these 
sub-systems in order to achieve effective thermal control. In many research and industrial 
applications it is necessary to perform complex thermal cycling. Such cycling may involve 
specific thermal ramp rates and extend over long periods of time (days) necessitating 
computer control. 

Even in a thermally controlled system where effective thermal insulation has been 
employed, effective thermal conduction to the sample has been used and sufficient heat is 
exchanged, this does not ensure thermal stability since typically the cooling system will not 
be continuously in operation and will have a certain time delay between activation and the 
onset of cooling. In practice therefore thermal control will consist of a feedback system of 
thermal sensors and thermal control which will introduce oscillation of the test section 
temperature. Additionally ensuring thermal stability may not ensure thermal uniformity as 
the application of cooling may not be physically at the same location as the source of heating 
or thermal loss which therefore leads to spatial temperature gradients. An obvious method 
to both stabilize the test section (or sample) temperature and achieve improved thermal 
uniformity is the use of a massive conductive (metal) test section element or sample 
mounting section. Here the thermal inertia of the mass achieves thermal stability and the 
high thermal conductivity of the mass ensures uniformity. It is however often difficult to 
find space to house such a massive element and the cooling/ heating time of this element 
will necessarily be long. In the case of the AWTSI and AWTSII facilities a compromise has 
been reached between the desired stability/ uniformity and the available space/ required 
response time. In the case of the AWTS-II facility at a temperature of around -120°C the 
uniformity of the test region is around 15°C/m and the stability is around ±2°C despite the 
use of a 0.1m thick sample (aluminum) mounting plate. 

8. Ice formation and sensing 

The transport of water vapor is often an important physical parameter for environmental 
simulators, both with regard to industrial and research applications. For example on Mars 
the desiccation of surface materials and subsequent frost formation (re-hydration) may lead 
to geophysical changes in the surface materials, specifically salt crust formation (which are 
widely observed) or even erosion. Similarly man made materials can be susceptible to 
weathering by the transport of water vapor from the surface. For the control of humidity at 
low temperatures it is necessary to both cool the sample and another element within the 
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environment. Typical research environmental chambers rely only on mounting the sample 
on a cold finger, in this case the sample will typically be at the lowest temperature within 
the environment and therefore attract out gassed water vapor and be hydrated (often 
forming surface ice). Other environmental chambers either cool the outer chamber or the 
atmosphere directly. In this case the sample will typically be warm compared to this 
environment and lose water vapor i.e. become desiccated. If the sample and environmental 
element(s) are thermally controlled independently then water diffusion to and from the 
sample can be regulated allowing detailed study of water transport phenomena. This does 
however require more sophisticated cooling and (computer) control systems. It is also then 
additionally desirable to have both cooling and heating systems for enhanced thermal 
regulation. Environmental (cooled) wind tunnels have been used occasionally in studies of 
snow transport, often involving processes of electrification which is relevant to the 
generation of electrical thunder storms and lightening [Maeno 1985, Schmidt et al. 1999]. 
This field of research is sparsely investigated and still largely empirical, despite being of 
great importance to human safety (anywhere that people interact with snow and ice) and 
relevant to understanding a variety of (dangerous) meteorological phenomena including the 
vast Arctic/ Antarctic regions. The technology for such research is readily available and 
straightforward to implement, involving the combination of a wind tunnel structure and 
(sub-zero) refrigeration techniques. 

9. Aeolian sand and dust transport 

Desertification involves the erosion and degradation of drylands which affects 25% of the 
Earths landmass and more than 2.1 billion people. Each year 12 million hectares of land are 
lost to desertification. This is the UN launched decade for deserts and the fight against 
desertification. The relevance of this work into Aeolian grain transport is to understand, 
predict and control this environmental effect [Ofori and Showstack 2010]. As well as being 
affected by the climate (and climate change) through changes in wind and surface 
conditions, Aeolian transport can also strongly affect the climate through the generation of 
aerosols (e.g. suspended dust). As an example, every year of the order of 1 billion tons of 
dust is removed from north Africa into the atmosphere and is carried great distances, 
around 200 million tons is subsequently deposited on south America. The affect of such 
aerosols in the atmosphere, also including smoke, clouds, ash etc., is complex. There is the 
need for a lot more research to understand the role of aerosols, including laboratory studies 
(environmental wind tunnel studies), remote and in-situ sensing as well as modeling. 
The definition of sand in the context of wind tunnel operation is granular material which 
cannot be suspended by the flow, but can be entrained i.e. can be removed from the surface 
by the wind shear. The definition of dust in this regard is therefore granular material which 
is fine enough to be suspended. Typically in nature sand particles will be transported by 
saltation wherein they are repetitively removed (entrained) from the surface, but cannot 
remain suspended and return to the surface in a ballistic impact referred to as a splash [Pye 
and Tsoar 1990]. This definition of sand and dust is not a strict size/ mass scale since it is 
dependent upon many physical parameters for example; the wind shear, atmospheric 
properties, gravitational conditions, etc. Under terrestrial conditions sand is typically 
greater than around 60|Lim up to mm size. Particulates larger/ more massive than this can be 
moved by wind shear without leaving the surface either by rolling, sliding or by impact of 
sand, this transport process is referred to as creep. 



Environmental Wind Tunnels 19 

Saltation is an effective erosion process due to the destructive nature of the splash. These 
impacts typically lead to chipping and the generation of dust sized particulates. Aeolian 
(wind driven) particulate transport is also effective at sorting granular material, specifically 
separating sand, dust and stones/ pebbles. This leads to the generation of the many sand 
dune forms observed, their character depending on the nature of the wind flow. The 
threshold wind conditions (surface wind shear) required for sand saltation transport has 
been widely studied, both in wind tunnel simulators and in nature following the pioneering 
work of Bagnold [Owen 1964, Greeley and Iversen 1985]. Despite this the work remains 
semi-empirical, mainly due to the complexity of the physical parameters involved, such as; 
wind induced lift and torque, adhesion, gravitation and the effect of splash. Similarly a large 
body of work exists in which the transport rate (above threshold) has been studied. Again 
resulting in an array of increasingly sophisticated and intricate semi-empirical expressions 
with especially the effect of splash being problematic [Shao and Raupach 1992, Rasmussen 
and Sorensen 2008]. Interesting aspects to be noted is that the rate of sand transport 
increases rapidly above the threshold wind speed (varying as the cube of the wind speed) 
and that due to the effect of splash saltation can be maintained below (by around 20%) the 
initial threshold speed once the process has begun. 

As discussed in the introduction the different gravity on other planets is problematic to 
reproduce in terrestrial laboratories. However, it is possible to vary the effective mass 
density either by using low density material or hollow structures (glass bubbles), this has 
been effectively utilized in wind detachment studies. Unfortunately a complex process such 
as the splash in saltation involves both the inertial mass (on impact) and the gravitational 
mass (during the trajectory) of the sand particulates. This makes this process impossible to 
entirely reproduce in a terrestrial laboratory and can probably only be studied through 
partial experimental simulation, modeling and observation in the respective gravitational 
field [Merrison et al. 2007]. 

Wind tunnels involved in sand transport are typically of the order of a square meter cross 
section, several meters long and are constructed as boundary layer simulators with the use 
of (upstream) flow control to reproduce the correct surface boundary layer (and shear stress) 
[Young 1989]. These flow control units include roughness arrays, inlet meshes, turbulence 
spires, and systems for injecting sand upstream (see figure 3). The results of laboratory 
based wind tunnels can be complemented by the use of field wind tunnels. These are 
portable wind tunnels which are taken into the field and thereby use real surface material 
and environments and possibly even the actual wind flow. 

Dust can be entrained and therefore transported by wind shear in two ways. The generally 
accepted process of dust entrainment on Earth is through the action of saltating sand grains 
impacting dust particulates and ejecting them into suspension. This process is extremely 
effective at wind speeds above the threshold for saltation. The entrainment of single (e.g. 
micron sized) dust grains directly by wind shear requires extremely high shear stress 
(extremely large wind speeds), however dust grains generally cohere and form larger 
aggregates, often up to mm sized. The mass density of these aggregates can be low (below 1 
g/cm3) and can therefore be mobilized at much lower wind stress than solid sand grains. 
Once detached from the surface the dust aggregates may disintegrate and liberate free dust 
grains into suspension. This process has been observed and quantified in (environmental) 
wind tunnel studies and helps to explain the otherwise paradoxical observation of dust 
transport in the Martian atmosphere (at wind stress below the sand transport threshold) 
[Merrison et al. 2007]. This transport mechanism may also be important for dust transport in 
areas not rich in sand, though dusty. 



20 Wind Tunnels 

An interesting erosion effect which has emerged as a result of studies in a Mars simulation 
wind tunnel is the possibility that saltation may also lead to alteration in mineralogy and 
thereby not just the generation of dust, but also for example oxidation. This may explain the 
presence of the reddish iron oxide giving Mars its distinctive hue. On earth mineral change 
due to presence of liquid water and high atmospheric oxygen content probably makes such 
erosion induced mineral change only a minor and not easily identified process, though may 
also have importance through the chemical reactivity of the erosion generated dust. The 
mineral change appears to occur through mechanical activation of freshly cleaved surfaces, 
for example in the case of silicate leading to an oxidizing surface which may be hazardous to 
organic material [Merrison et al. 2010]. This phenomenon is as yet poorly researched and 
requires the further application of laboratory simulation. 

10. Aerosol formation and sensing 

Open circuit wind tunnels are poorly suited to studies of suspended particulates, i.e. 
aerosols, due to problems of contaminating the environment with the aerosol particulates 
and also since the aerosols only perform a single pass through the detection volume. This 
limits the type of aerosol dynamics which can be studied, for example investigating their 
temporal evolution would not be possible. In a re-circulating wind tunnel the aerosol can be 
confined and studied for long periods of time, limited only by loss due to settling and 
adhesion [Merrison et al. 2002, Merrison et al. 2008]. 

As discussed previously fine solid particulates (i.e. dust grains) generally tend to cohere and 
are found as large aggregates. In order to inject dust particulates into suspension these 
aggregates must be dispersed. A problem affecting many dust injection systems 
(aerosolizers) is that of blockages caused by dust aggregation (clumping). A system which 
functions well and is widely used in the medical industry for dry powder aerosol 
dispersion, is the use of a gas jet. When merged with the dust material, the process of gas 
expansion (usually from a nozzle) disperses and separates the aggregates. This same 
technique has been successfully employed in planetary aerosol simulation where (relatively) 
high pressure gas is passed through a chamber containing the dust material and 
subsequently expands at high velocity into the wind tunnel flow. This system of gas 
injection is also successfully applied in the case of liquid (droplet) aerosols. Liquid aerosols 
are commonly used in pharmaceutical applications, for example inhalers and the many 
types of commercial 'aerosol cans'. Another method for entrainment (suspension) which 
resembles that seen in nature, involves emplacing a thick dust layer in the wind tunnel and 
increasing the wind speed briefly in order to suspend it. Aerosols may also be formed as a 
result of condensation/ precipitation from the gas phase. In this case a system of vapor 
injection from an external source may be utilized. Alternatively the vapor source could be 
internal to the system and be generated by heating to generate sublimation/ boiling. This 
may be especially effective if the process of condensation (or nucleation) is to be studied. 
Surprisingly there is still a great deal yet to be understood about the properties of ice/ water 
aerosols (e.g. clouds and snow). For example; the electrification process or processes which 
are responsible for generating lightening, the details of nucleation (for example by dust, 
radiation or ions), the micro/ nano structure and adhesive/ cohesive properties. It seems 
likely that this will be an active research field in the near future and that this will involve the 
construction of new low temperature environmental simulators, possible including re- 
circulating wind tunnel systems. 
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Particulate deposition from suspension is of great importance with regard to the climate as 
well as being a hazard to instrumentation, machinery and human health. Lacking commercial 
instrumentation to quantify this deposition process, a prototype optoelectronic instrument has 
been developed to determine the deposition rate of particulates [Merrison 2006]. This 
instrument also employs electrodes to generate electric fields and thereby attract electrified 
aerosol particulates. This method of determining the electrical charge (and sign) of suspended 
particulates could be of great use in industries dealing with aerosols. Electrostatic fields are 
already used routinely in industrial applications to remove or control particulates both in 
suspension and otherwise. Granular electrification is a ubiquitous effect and must be 
considered wherever granular material is found [Rasmussen et al. 2009]. In the case of dust 
particles electrification leads to the formation of aggregates which is key to the process of 
entrainment at low wind speeds (in the absence of saltating sand) [Merrison 2004a]. 
On earth aerosols may include (or consist) of biological material or micro organisms (such as 
pollen, virus, bacteria, spores, etc.). Bio-aerosols of this kind are clearly of great 
environmental interest with regard to health and ecology. It would be of much interest to 
study bio-aerosols and their production/ transport mechanisms in a controlled 
environmental wind tunnel. Although such research is not currently active, there are plans 
to begin studies of this kind with relation to planetary protection (i.e. the protection of other 
planets from terrestrial micro organisms) by the space agencies NASA and ESA. 
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1. Introduction 

Investigation of bodies liable to high speed or variable airflows is a subject of importance. 
Determination of force and torque vectors acting on the latter is necessary to establish 
strength of the body and flow behavior. Airflow influence on a body (or system of bodies) is 
used in science and industry to define the values of air resistance parameters acting on this 
body. Resisting forces and torques, lifting or carrying forces and torques, and air attacking 
angles are examples of such types of parameters. This work deals 6-DOF force and torque 
measuring suited for use in air tunnels, comparing it with the existing methods used for the 
same purpose. Except for the body under investigation and its fastening elements nothing 
else is placed inside the tunnel. Six force sensors (load cells) and their fastening design is 
organized outside the air tunnel by means of a metal cubic (or other) frame were suspended 
on a massive base. Mathematical expressions describing the measurement and computation 
process of the values of the above mentioned 6 force and torque vectors are proposed. 
Examples of measuring processes and graphical presentations of some experimental results 
are given in this work. In this work we propose further development of the discussed design 
by reducing the mass and inertia moments of this frame and relative positioning the model 
to it, in such a way reducing the caused by these dynamic parameters disturbing the 
measurements. In this case the model pattern for air flow influence measurements were 
made on a flat disc-like body. 

2. Brief overview of wind tunnels. 

Means for studying and determining the forces and torques acting on bodies influenced by 
high air flow speeds or subjected to variable wind speeds in air tunnels are required for both 
science and industry. This airflows influence on a body (or a system of bodies) is used for 
estimation the air resistance acting on these bodies or body. In the authors work [1] is 
brought a brief list of three balances that are usually applied for measuring wind drag forces 
and torques, lifting or carrying forces and torques, "air attacking angle" are examples of such 
parameters. 

Partly the described in the mentioned paper measuring lay-outs are illustrated by examples 
given in Ewalds paper [2] from 2000. Another type of internal six- component wind-tunnel 
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balance is that having three-dimensional mechanism, as described by Gorlin and Slezinger 
[3], 2000. A solutin alike to mentioned here belongs also to Corliss and Cole [4] 1998. In these 
cases a multi-component balances are placed in the tunnel being by means of elastic cables 
connected to corresponding elements of the measured pattern. This set disturbs the air flow 
in the tunnel to some extent. 

The model can be via an enough elastic rod connected to the tail part of the investigated 
model and provided with differently oriented strain gauges gives the needed measured 
information describing the acting forces and torques. Pope and Harper [5] 1999 describe 
developed some calibration rigs fir internal strain gauge balances. 

The closest to the proposed design is the model placed on a platform inside the tunnel. And 
the forces are measured with 6-DOF hexapod-like balances placed outside the tunnel. (See 
Nguen et al.[6], 1002 and Kerr et al [7],1989) 

In this work, we propose to use a fourth orthogonal structure of a parallel robot directly in 
the chain of force and moment components, acting on the model, and computation. An idea 
came to shorten the process of investigation of the forces and moments acting on the model 
in the tunnel, which avoids the amount of unnecessary bodies in the tunnel and their 
interference. This proposed measurement method in the wind tunnel is close to the discussed 
earlier external 6-DOF to certain extent improved approach. Some ideas of this kind where 
discussed in the authors work Chapsky et al. [10], 2007, and the work [11] of Lui, S.A. and 
Tso, H. L., 2002. The main purpose of this paper is to propose and demonstrate 
experimentation with a wind tunnel structure including the model to be investigated, a 
measuring element combined with a computer. This system (shown in Fig.l.) uses the 
corresponding functions (or equations) for each of six force and moment components 
depending upon the air speed and the shape of the investigated model. This work is 
devoted to introduce a different measuring approach and shows its advantages offering the 
fourth kind of dealing with the problem. The fourth kind of balance proposed in this paper 
has once calibrated load cells and once written computation program for any fixed model 
used for investigation. 

only one rod is placed in the tunnel for fastening the investigation pattern; 

the balances and measuring sensors (load-cells) are located outside the tunnel 

supplying the needed information to the computer; 

values of the acting, due to the wind, forces and torques are automatically shown on the 

computer's screen; 

there are no interactions between the six measured and computed values as it often 

happens in other systems. 
Later we show 7 pages with computed examples of all here mentioned types (even those 
which we are not able to realize for our tunnel). 

3. Experimental arrangement 

A general view of the wind tunnel built at the Mechanical Engineering Department of Ben 
Gurion University of the Negev, Israel is shown in Fig.l. 

In Fig.2 a close view of the wind tunnel near the outer frames and the places for movable 
frame's suspension are given. Also the door for the pattern fastening and its orientation is 
seen here. Fig.3 shows the load-cells used in this design. It used for the measuring the value 
of the mechanical response to the wind forces. The maximum air speed in the tunnel is 40 
m/sec. Length of the tunnel is about 30 m. 
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Fig. 1. Wind tunnel device in the aero-dynamic laboratory of BGU of the Negev, (general 
view). 




Movable frame 



Immovable 
frame 



Fig. 2. Wind tunnel. The movable aluminum frame provides the information about the X, Y, 
and Z displacements of the pattern and its rotations around these axes. The steal made 
immovable frame is the base of the arrangement. 
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(a) (b) 

Fig. 3. The Data-Logger of the 34970 type: (a) S-beam load cell, (b) view of the suspension 
including the Data-Logger. 

Cross-section of the tunnel equals 0.7 x 0.7 m^. The pattern is fixed in the middle of the 
earlier mentioned rod as it is schematically shown in Fig.4. This experimental case discussed 
in work [1]; 




Fig. 4. Lay-out of the movable old frame primarily designed and built around the tunnel: a = 
470, b = 1050, c = 1030, L = 200, (the dimensions are millimetres). 
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Fig. 5. Pattern location in the wind tunnel: flat disc placed in it. 




Fig. 6. Relative design of the new and old (previous) frame construction. 



4. Main idea of the new frame design. 

To improve the dynamical performance of the proposed apparatus the set of different frame 
arrangements has been investigated. The difference between experimental arrangements lies 
in the mass of the entire design, in the position of the tested pattern (flat disc) relative to the 
frame, and in the symmetry of the frame shape. 

To reduce the mass of the frame without changing the geometric dimensions the new 
arrangement has been proposed. Its design is clear from figure 6. In the Fig.6 the red frame 
is the new proposed design, the "gray" profile 5 demonstrates the wind tunnel, the 'black" 
frame 4 is the previous design, 1, 2, 3 are the fastening points of the frame suspensions for 
both cases (new and old design), the green sphere 6 is the investigated pattern. 
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Every point of the frame suspensions has two directions of force measuring. For instance, 

the point 1 measures in directions "x" and "z", the point 2 measures in directions "z" and "y", 

point 3 in directions "x" and "y". 

The investigated models were conventionally named as follows: 

"Old'' - with full set of 12 edges, ''New'' - with reduced mass (4 edges); 

"Regular" - with pattern, located outside the frame; 

"Central" - with pattern in the geometrical center of the frame; 

"Symmetric" - with cube-shaped frame coverage figure. 

The investigated frame designs and their dimensions are shown in Table 1. 



Old Regular 



New Regular 



Old Central 



New Central 



Old symmetric 



New Symmetric 







-<u. 





a = 4:70 
b = 1050 
c = 1030 
L = 200 



a = 1050 
b = 1050 
c = 1050 
L = 200 



a = 470 
b = 1050 
c = 1030 
L = -a/2 



a = 470 
b = 1050 
c = 1030 
L = -a/2 



a = 1050 
b = 1050 
c = 1050 
L = -a/2 



a = 1050 
b = 1050 
c = 1050 
L = -a/2 



Table 1. Proposed frames, their design and dimensions. 

All these "Old" and "New" kinds of frames are checked dynamically and therefore can be 
accordingly analyzed. 

The item we propose for publishing is somehow continuation of the paper published earlier 
and caused your attention. We speak about an improved design of the same principal 
structure, however improved because its main mass is considerably smaller despite the 
same geometric dimensions. The said is clear from the figure 6. 



5. Dynamical investigation 

Comparative dynamic investigations of virtual models of tunnel balances have been carried 

out by means of the interactive computer-based motion simulation software MSC 

ADAMS/View. 

Every virtual model was loaded consequentially in all directions by harmonic forces or 

moments U(t) with varying frequency. 

U(t) = Uo sin (g) (t)*t) = Uo sin ((1000*t)*t), 

Where Uo - initial magnitude of the force or moment, oo (t) -frequency, which depends on 
the time according equation: oo (t) = 1000*time. 

The responses of each spring have been measured and Bode plots were calculated with 
purpose to define the bandwidth of tested model. Bandwidth corresponds to the 
frequencies, where Bode-magnitude plot lay in the ±3dB range relative to the initial (steady- 
state) value. For convenience, the Bode plots were shifted to zero value in initial state by 
multiplying the measured results on some factors. All factors were calculated by symbolic 
mathematics software such as Mathematica-7. The applied Mathematica-7 program is 
shown in appendix 1. 
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Old Regular (Py, Spring 2) 
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Frequency (Hz) 
End Time 1sec, Steps 4000 Py = Sin((1000*time)*time) 



a 



1 1 1 mill iiiiiini II 



1 



10000.0 Q-l 1-0 yiO.O 100.0 lOOO.OIOOOO.O] 

4.88 Frequency (Hi) 




4.881/sec 



New Regular (My, Spring 3) 
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59.57.281/sec 



New Center (My,Spring 1) 



100.0 1000.0 10000.0 



Frequency (Hz) 
, Steps 4000 My=Sin((1000'time)'time) 




1.0 10.0 lOOlK^OOO.0 10000.0 
Frequency (Hz) 




137.701/ sec 



Old Symmetry My, Spring 2) 



g-20.0' 



a 



-<: 



10.0 100.0 1000.0 10000.0 



Frequency (Hz) 
End Time Isec, Steps 4000 My = Sin((1000'time)*time) 



0.1 1.0 10.0 IOO.CM^O.O 10000.0 

Frequency (Hz) ,^ 




108.401/sec 



II _, 



New Symmetry (Px, Spring 4) 

fz-1.5 



10.0 100.0 1000.0 10000 

Frequency (Hz) 



1.0 10.0 100.0^00.0 10000.0 
Frequency (Hz) 215.82 




End Time Isec, Steps 4000 Px = Sin((100D'time)'time) 



215.821/ sec 



Table 2. The determinative Bandwidths 
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The worst Bode plots relative to spring forces and corresponding bandwidths for every 
force and torque directions (Px, Py, Pz, Mx, My, Mx, according notations in Figure 5) are 
shown in Appendix 2. The determinative ones for each design are shown in Table 2. 
As clear from Table 2, the most desirable is the design ''New Symmetry'', which has 
bandwidth 215.82 Hz. The "New central" is applicable too. If it is not possible to locate the 
pattern in the center of the frame due to the limitations of the wind tunnel design the "New 
regular" model can be applied, but with bandwidth 113.28 Hz. All these models are better 
than earlier model "Old Regular", which has bandwidth 4.88Hz in Py and Pz directions, 
despite the fact that in Px direction it has bandwidth 156.25Hz. 

6. Conclusions 

1. The proposed device requires a universal sequence and simple action for receiving the 
measured data. Investigator must only fasten the pattern on the mentioned rod. 

2. Dynamic estimations of the virtual model of proposed tunnel balance by means of the 
interactive computer-based ADAMS/ View system showed that the working frequency 
of the proposed method and device are limited to a -100^200 Hz bandwidth. 

3. The translational isotropy of the proposed device is defined by the independence of the 
sensitivity of measurement from the direction of the operating forces. 

4. For estimation of the anisotropy of devices, the anisotropy index which equals to the 
ratio of maximum and minimum stiffness values is applied. Two types of indexes are 
separately used: for translational and rotational stiffness values (see Table 2). 

5. The proposed devices "New Symmetry" and "Old Symmetry" has better dynamic features 
than conventionally used systems: it is fully isotropic from the point of view of translational 
stiffness and has a high level of isotropy from the view point of rotational stiffness. 

6. The proposed device shows at least the theoretical possibility to improve the dynamic 
properties of the wind tunnel comparing with the already offered. 
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Appendix 1. Mathematica-7 computation program. 

^Removing all data"- 
Remove ["Global +"] ; 




^I^ Old Fram€ 



"Set of equations of static equilibrium of forces and moments' 

Ml = {Px+fT4 +1^5 ==0, Py +n2 +£ie » 0, Pz +^1 +£I3 ==0, 

Mx+N3*b/2-Hl*b/2-N2*c/2+n6*c/2==0, 

My+N3*L+Nl+(L + a) -N5+c/2+H4*c/2 ==0, 

Mz+N5*b/2-H4*b/2-N6+L-lI2*(L + a) ==0}; 
MatrixForm [%] 



ET4 +N5 +Px == 
112 + K6 + Py == 
HI +K3 +Pz ==0 

+ + + t ] K == 

2 2 2 2 

i£i . ^ +H3 L+Hl <^ +L) +My == 
2 2 

.^ + ^.N6L-N2 (a+L) +Mz==0 
2 2 
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"Solving the system of equations" 

Sol vMatrix = Solve [Ml, {HI, W2 , H3 , N4 , H5 , Ne>] 

2aMx-2bMy-2cME-acPy-2cLPY+3abPz+2bLPz 



{{n3. 



Nl -^ -■ 



H4 



N5 



N2 



4 ab 
-2aMK + 2bMY+2cMz+acPy+2cLPy + abPs-2bLPz 

4 ab 
2aMx+2bMy-2GM3+2bePx-acPy-2c!LPy-abP3-2bLP2 

4 b c 
-2atlx-2bHy+2cMs-H2bcPx+acPy + 2cLPy + abPs+2bLPa 

4 b e 
-2aMx-2bMy-2cMz+acPy-2cLPy + abPz+2bLPz 



N6 -* -■ 



4 a c 
2aMx-H2bMy + 2cHz+3acPy+2cLPy-abPz-2bLPz 



4 ac 



}} 



"Conditions of the simulation" ; 



condPx = {Py ->0, Pz-iO, Mx-^0, My-^0, Mz->0} 
condPy = {Px ->0, Pz-»0, Mx-^0, My-^0, Mz-^0} 
c 
c 



=ondPz = <Py ->0, Px-*0, Mx->0, My-»0, Mz-*0> 
riondHx = {Py ->0, Pz-»0, Px-»0, Hy-^0, Mz-*0} 
condMy = {Py ->0, Pz-»0, Ms-»0, Px-»0, Mz-»0} 
condMz = {Py -> , Pz -> , Mx -* , My -» , Px -» 0} 



"Geometric conditions: Regular,, Central, Symmetric, 
and L for central and symmetric configurations"; 

condGeomReg = {a -> 470 , b -> 1050 ,, c -» 1030 , L -> 200} ; 
condGeomSym = {a -> 1050, b -> 1050, c -* 1050); 
condGeomCentr = {a -> 470, h -> 1050, c -» 1030}; 
condL = {L -* -a /2}; 
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"Numerical results"; 

R»gPx ■ SolvMwtriK / , conclPx / , condGeomReg 

!(ti3 ^ , HI ^ , H4 ^ , N5 ^ , N2 ^ , 116 ^ oU 

R»gPy ■ SoivHfltriX / , condPy / . condGeomReg 

ff 2987 Py 2987 Py 29 Py 29 Py 7 Py 181 Py^^ 

J {n3 ^ , HI ^ , N4 ^ , M5 ^ , N2 ^ , H6 ^ \ \ 

Li 5580 6580 140 140 188 188 ^^ 

L P2 

RegP; ■ SolvMatriK / . condPz / . coridG«oinR«g 

.- 181 Pi 7Pz 87 Pz 87 P= 9135 Pi 9135Pz.. 

\\n3 ^ , HI ^ , H4 ^ , H5 ^ , N2 ^ ■ , N6 ^ \\^ LPs 

'•'• 188 188 412 412 19364 19 364 ^ J 



Rft^x o SolvMatrix / . condMx / . condGeomReg 

«Mx Mx 47 Mx 47 Mx Mx Mx . . 
N3 ^ , HI -^ , N4 ^ , H5 ^ , H2 -* , N6 ^ [ [ 
2X00 2100 216300 216300 2060 2060'''' 

Re^y = SolvMatrix / . eondMy / . oondGeomReg 

«My My My My 21 My 21 My.. 
N3 ^ , HI -^ ^ , N4 -» , H5 -* , H2 ^ , H6 ^ [ \ 
940 940 2060 2060 19364 19364''-' 

Re^z = SalvMatrix / . condMz / . condGeomReg 

«103Mz 103 Mz Mz Mz Mz Mz . . 
N3 ^ , HI -* , H4 ^ , H5 ^ , H2 ^ , H6 -» ^ [ \ 
9S700 98700 2100 2100 940 940-''' 

CetitrPx = SolvMatrix / . condPx / . oondL / . condGeomCentr 

||h3 ^ , HI -> , M4 -^ , N5 -* , H2 -» , H6 ^ o| j 

CentrPy = SolvMatrix / . condPy / . condL / . condGeomCentr 
H3^0, Hl^O, M4^0, H5^0, H2-» , H6-* \\ 

CentrPz = SolvMatrix /. condPz /. condL /. condGeomCentr 

rr ^^ ^= n 

||tl3 ^ , HI -> , H4 -* , N5 -* , H2 ^ , 116 ^ j| 

CentrMx = SolvMatrix / , condMx / . condL / . condGeomCentr 

„ Mx Mx 47 Mx 47 Mx Mx Mx . . 

J Jtl3 -* , Nl -* , H4 -* , H5 ^ , N2 -^ , N6 -* [} 

'■'■ 2100 2100 216300 216300 2060 2060-'-' 

CentrMy = SolvMatrix / , condMy / . condL / . condGeomCentr 

«My My My My 21 My 21 My.. 
H3 -* , HI -» , N4 -» , H5 -^ , H2 -> , N6 -» \ \ 
940 940 2060 2060 19364 19364-'-' 

CentrMz = SolvMatrix / . condMz / . condL / . condGeomCentr 

«103Mz 103 Mz Mz Mz Mz Mz . . 
H3 -* , HI -^ , H4 -* , H5 -4 , H2 -» , H6 -» \ \ 
98 700 98 700 2100 2100 940 940-'-' 
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SymPx = SolvMatrix /. condPx /. condL j . condGeomSYm 

|[n3 -» , Nl -* , N4 -• ■, N5 -. ^ , N2 -* , N6 ^ o|| 

SymPy = SolvMatrix / . condPy / . condL / . condGeomSym 
N3 -. , Nl -I , N4 -• , N5 -. , N2 -• , N6 -• 1 j 

SymP z = SolvMatrix /. condPz /. condL /. condGeomSym 

UPz Pz .. 

N3 -. -, Nl -• , N4 -• , H5 -» , N2 -• , N6 ^ 0| j 

SyniMx = SolvMatrix / . condMx / . condL / . condGeomSym 

UMx Mx Mx Mx Mx Mx ^ ^ 
N3 -. , Nl -• , N4 -• , H5 -• , N2 -• , N6 -• [ 
2100 2100 2100 2100 2100 2100'''' 

SymMy a SolvMatrix /. condMy /. oondL / . condGeomSym 

UMy My My My My My ^ ^ 
N3 -» , Nl -• , N4 -• , H5 -. , N2 -• , N6 -. f 
2100 2100 2100 2100 2100 2100''-' 

SymMz = SolvMatrix / . condMz / . condL / . condGeomSym 

UMz Mz Mz Mz Mz ^^ ^^ 
N3 -. , Nl -• , N4 -• , N5 -• , N2 -• , N6 -i f [ 
2100 2100 2100 2100 2100 2100''-' 

RegPxNevr s SolvMatrix /. condPx / . condGeomSym /.!,-* 200 

Px Px 

— ^ H5 -. 

2 2 



{[n3 -.0, Nl-iO, N4-. , N5-. , N2-i0, N6 ^ o}| 



RegPyNew = SolvMatrix /. condPy / . condGeomSym / . L -i 200 

U29 Py 29 Py 29 Py 29 Py 13 Py 71 Py.. 
N3 -. , Nl -» , n4 ^ , W5 ^ , N2 ^ , N6 -> \ 
84 84 84 84 84 84 '''' 

RegPzNew = SolvMatrix /. condP z /. condGeomSym / . L -i 200 

U71PZ 13Pz 29Pz 29Pz 29 Pz 29Pz.. 
N3 -. , Nl -• , N4 -• , N5 -• , N2 -» , N6 ^ [ 
84 84 84 84 84 84 -■ -" 

Re^xNew = SolvMatrix /. condMx /. condGeomSym / . L hi 200 

UMx Mx Mx Mx Mx Mx ^ ^ 
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Appendix 2. Bode plots of the forces and torques acts on the pattern relative 
to spring forces. 
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1. Introduction 

As is well known, a wire-driven parallel manipulator is a manipulator whose end-effector is 
driven by a number of cables instead of rigid links. It shows several promising advantages 
over its rigid-link counterpart, such as simple light-weight mechanical structure, low 
moment inertia, large reachable workspace and high-speed motion. In the 1980s, the 
National Institute of Standards and Technology (NIST) in America invented a wire-driven 
parallel manipulator named RoboCrane for shipyards (Albus et al, 1993). So far, wire-driven 
parallel manipulators have been applied in load lifting, industrial machining, virtual reality 
and astronomic observation (Dekker et al, 2006; Ning et al, 2006; Ma & Diao, 2005). Because 
of the advantages and unique features of wires, wire-driven parallel manipulators have 
attracted a great attention in robotics literature. The first general classification was given by 
Ming and Higuchi (Ming and Higuchi, 1994). Based on the number of wires (m) and the 
number of degrees of freedom (n), wire-driven parallel manipulators were classified into 
three categories, i.e. the incompletely restrained positioning mechanisms {m<n+l), the 
completely restrained positioning mechanisms {m=n+l) and the redundantly restrained 
positioning mechanisms (m>n+l). Yamamoto et al. presented basic dynamics equations and 
a feedback control method based on exact linearization for the incompletely restrained 
positioning mechanisms (Yamamoto et al, 2004). Hithoshi et al. studied a robust PD control 
using adaptive compensation for translational wire-driven parallel manipulators of a 
completely restrained type (Hithoshi et al, 2007). Zi Bin et al. developed a fuzzy plus 
proportional-integral control method for the cable-cabin mechanism of 500m aperture 
spherical radio telescope (Zi et al, 2008). Yu Kun considered active stiffness control schemes 
as optimization problem with different criteria for redundantly restrained positioning 
mechanisms (Yu, 2008). In essence, a wire-driven parallel manipulator can be considered as 
a complex, time-varying, strong-coupled, multiple input and multiple output, and nonlinear 
system. Since the wires can only pull and not push on the platform, dynamics and control 
are key issues for high-precision motion of wire-driven parallel manipulators. 
Wind tunnel tests of aircraft models are widely utilized to investigate the potential flight 
dynamics and aerodynamic characteristics of aircrafts at their early developing stage. Wire- 
driven parallel manipulators have been introduced to wind tunnels as flexible suspension 
systems of aircraft models in recent years (Liu et al, 2004). The posture of the scale model 
corresponding to the stream line of airflows can be adjusted by controlling the length of 
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wires to implement the six degree-of-freedom free flight motion. The aerodynamic forces 
exerted on the scale model can be calculated by measuring the tension of each wire. 
Comparing with traditional frame suspension systems, wire-driven parallel manipulators 
for wind tunnels have advantages in less aerodynamic interference and high precision of the 
test results. Preliminary achievements have been made in the Suspension ACtive pour 
Soufflerie (SACSO) project about the wire-driven parallel suspension system in low wind 
tunnels sponsored by Office National d'Etudes et de Recherches Aerospatiales (ONERA). 
The achievements include architecture design, workspace computation, force control and 
build-up of a prototype of the wire-driven parallel manipulator (Lafourcade, 2004). Zheng 
Yaqing et al. have developed some fundamental theoretical research work on workspace, 
wire tension distribution, stiffness, kinematics and control of the manipulators. Because of 
weak stiffness of wires, the aircraft model would deviate from the planned trajectory when 
it is in the streamline flow. The trajectory errors have significant effect on the force and 
moment measurement. Hence one challenging issue is to accurately implement the attitude 
control for wire-driven parallel manipulators in wind tunnels. 

The flexible suspension system in wind tunnels proposed by Zheng, which can be viewed as 
a six degree-of-freedom eight wires driven parallel manipulator, is investigated in this paper 
(Zheng, 2004). In order to decrease the trajectory errors and improve the measurement 
precision, it is necessary to enhance stiffness of the flexible suspension system. In case of 
wire-driven parallel manipulators with redundant actuations, the stiffness of the 
manipulators have been researched by Yu (Yu, 2008) and Saeed Behzadipou (Behzadipour & 
Khajepour, 2006) respectively, based on the stiffness definition and the equivalent spring 
model. In this paper, an analytic expression of the stiffness of the flexible suspension system 
in wind tunnels is derived by using the differential transformation principle. In order to 
hurdle a low rigidity and poor positioning accuracy caused by the minimum wire tension 
solution, an optimal tension distribution method is applied for the enhancement of stiffness 
in lift, along- wind and pitching directions. The method resolves the uncertainty of wire 
tensions of the suspension system. 

The motion control of the flexible suspension system in wind tunnels can be realized either 
in end-effector space or in joint space. The pose of the aircraft model must be measured in 
real time during the former control process. Measuring the pose of the aircraft model in 
wind tunnels is rather challenging, because the cross section of wind tunnels is limited and 
the existence of equipments disturbs air flows. Moreover, it is not desirable to obtain the 
pose of the aircraft model using direct kinematics, because of lots of time required by 
complicated calculation. Hence, a computed torque controller in joint space is employed for 
the flexible suspension system in wind tunnels. A dynamics compensation is introduced to a 
conventional proportional differential controller, so a modified proportional differential 
control strategy in the wire length coordinates is developed based on stiffness enhancement. 

2. System description 

Figure 1 shows the flexible suspension system driven by eight wires. Each wire is attached 
to the aircraft model at one end, and threads the pulleys mounted to the wind tunnel and 
winds around an actuated reel at the other end. The actuated reels allow the control of the 
pose of the aircraft model by controlling the length of their respective wires. The 
aerodynamic loads on the aircraft model can be calculated through measuring the wire 
tension by strain gages. 
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Fig. 1. The flexible suspension system for wind tunnel 
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Fig. 2. Geometric definition of the suspension system 

All geometric quantities are shown in Fig. 2. OXYZ and O'X'Y'Z' are coordinate frames 
attached to the wind tunnel and the aircraft model, respectively. C is the center of gravity of 
the aircraft model. The point where the z'th wire leaves the reel is denoted by Bi, and the 
connecting point on the aircraft model is denoted by Pi. The rotation matrix of the O'X'Y'Z' 
with respect to OXYZ is represented by 



°R^, = rot{z,y/)rot{y,^)rot(x,^) 

cos^cosy^ COS ^ sin y^ sin ^- sin ^ COS ^ cos ^ sin y^ cos ^ + sin ^ sin ^ 

sin^cosy^ sin ^ sin y^ sin ^ + cos ^ cos ^ sin ^ sin y^ cos ^- cos ^ sin ^ 

-siny^ cos y^ sin ^ cos y^ cos ^ 



(1) 



where ^ , /3 and I// are the roll, pitch and yaw angles of the aircraft model respectively. 
The length of the /* wire is expressed by 



h =||%|L = VC^/ - "Po' - "Ro'^'PifCB, - °Po' - "Ro^'P,) for;=l,2,...,8 



(2) 
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where 'L^ = "B^ - "P^, - 'R^, "'/^. , 

^P^, =[^x^, "y^, ^z^,Y is the position vector of the mobile frame's origin, 

"" /^ = [""'xp "' yp "" z^ ]^ is the position vector of point Fi in the mobile frame O'X'Y'Z', 

°B- =[^Xg "yg ^Zg ]^ is the position vector of point Bi in the fixed frame OXYZ. 

Differentiating Eq.(2) with respect to time, and then assembling the eight resulting equations 
into matrix form, we obtain 



/ = -rx 

where I = [li I2 '" /g]^ is the wire-length vector, 

^ = ["^o' "yo' "^o' ^ P W^ is the posture vector of the aircraft model. 



(3) 



^U: = ^L- Yh\ is the unit vector alone the z't^ wire. 

I 1 1 II ;||2 o 

The equation of static equilibrium can be written as 

JT + F = 

where T = {t^t2 -" t^ Y is the wire tension vector, F -■ 
torques acting on the aircraft model. 



: R^""^ is a pose-dependent matrix. 



M, 



(4) 



summarizes all other force and 



3. Analytic stiffness 

The influence of the wire tension on stiffness of the flexible suspension system is 
investigated, and an analytic expression of the stiffness is derived from the differential 
transformation principle. When an infinitesimal wrench dF is applied to the aircraft model, 
the posture of the aircraft model changes by an infinitesimal deflection dX . The Stiffness 
matrix K of the suspension system is 



K 



' dx' 



dX dX 



(5) 



For the first term in the equation (5), dJ can be expressed by the product of an infinitesimal 
deflection dX and a three-dimensional matrix which excludes dX . Assuming the matrix H 

is equal to — , we obtain 

^ dX 



H = [H^ H^ 



H,]eR' 



H, 



-I 'p. X 

-CRy'P,)x rL.x][CRy'P.)x]-[CRy'p.)xrP.x]_ 



R' 



(6) 



where ( ) x is the operator representing cross product. 
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As for the second term in the equation (5), we have 

-/— = -j^.^ = kJdiag{—)J^ = kJdiag{lT\\ + k-\))J^ for z=l,. . .,8 (7) 

dX dL dX Iq- 

where k = EA , 

E is Young's modulus of a wire, 

A is the cross section area of a wire, 

/^ is the currently measured length of the z'th wire, 

/q^ is the original length of the z'th wire. 

It is pointed out that the contribution of ^ ~ t^ to the stiffness of the suspension system can 

be neglected because it is much less than one. The stiffness of the suspension system consists 
of two parts, while the first one is mainly influenced by the wire tension and the other one 
depends on geometrical arrangement of the wires and posture of the aircraft model. 
Supposing the external wrench F acted on the aircraft model is known, the wire tension in 
equation (4) can be written as 

T = -JF + Null(J)X (8) 

where J^ =J^(JJ^y^ eR^""^ is the Moore-Penrose inverse of matrix /, Null(J)eR^''^ is a 
matrix whose columns form a basis for the null-space of matrix /, A = [Xy A2Y eR^""^ is a 
column vector of two arbitrary real numbers. 

The solution in equation (8) consists of two parts: the first one is the term -J^F , which 
represents the minimum-norm solution that minimizes the 2-norm ||r|| . The second part 
Null{J)A, is an arbitrary vector in the mull-space of matrix / and, affects the distribution of 
the wire tension without affecting the force and moment at the aircraft model. Equation (5) 
can be rewritten as 

K « H{JF - Null{J)X) + k'Jdiag{-)J^ (9) 

It is clearly seen that the wire tension can be changed by adjusting the two elements of the 
column vector X , and then the wire tension can make an impact on the stiffness of the 
system. 

4. Dynamic models 

4.1 Dynamic Model of the aircraft model 

By using Newton-Euler's laws, the motion equations of the aircraft model can be written in 
the following form 

8 

mx + mco X ""C + mco x(o)x ""€) = ^ ""u-t- +mg + F^ 

(10) 

i''Cxx + Ia) + m{(0 X ^C) X X + 6? X {Im) = Yj"Pi^ "^th + ''Cxmg + M^ 
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where x = ["x^,"y^,"z^,f' represents the linear velocity of the reference point O' of the 

aircraft model, 

co = [^ fi y/Y is the angular velocity of the aircraft model, 

m is the mass of the aircraft model, 

^ = [00 gY and scalar g is the gravity acceleration, 

^C = ^R^< ^ C and ^ C is position vector of the center of gravity of the aircraft model in the 

mobile frame O'X'Y'Z' 



I=''RJ/Rl and /„, 



is the inertia tensor of the aircraft model in 



the mobile frame O'X'Y'Z', 

F^ and M^ are the force and moment exerted by aerodynamic load on the aircraft model. 

Equation (10) can be re-written into a compact form as 



where m 



N = 



W^ 



w^ = 



M{X)X + N{X, X)X = W^+Wg+JT (11) 

and / e R^^^ is the identity matrix, 

and G R^""^ is the zero matrix, 

is the aerodynamic wrench acted on the aircraft model, 

is the gravity wrench exerted on the reference point O' of the aircraft model, 

is the velocity vector of the aircraft model. 



ml -m^Cx 
m^'Cx I 

-m((Ox''C) 

m(cox^C)x -Id) 

Me 

mg 

""Cxmg 

X 
CD 



4.2 Dynamic model of the drive units 

A drive unit is composed of a motor, a gear reducer and a winch. The dynamic equation of 
the drive units is given as follows 



AO + Ce + rT = T 



(12) 



with A = diag{a]^,a2,---ag) , a-=a^-+— ^, 

n 



C = diag(c^,C2,'--c^),c.=c^ 

T 
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where a^ ■ , c^^ denote the moment of inertia and the vicious friction coefficient of the z'th 

motor, 

a^^ , c^i denote the moment of inertia and the vicious friction coefficient of the ft^ reducer 

and winch, 

r^^ is the radius of the z'th winch, 

n is the reduced ratio of each gear reducer, 

0- is the rotational angle of the z't^ motor, 

r- is the output torque of the z't^ motor. 

4.3 The elastic model of the wires 

The relationship between the change of the wire length and the rotational angles of the 
motors is 



= r-U 



'B, 



""B, 



'B, 



'B, 



'r.'pA 



'P,-'R,'PJ 



(13) 



Successive time derivatives of equation (13) yield 



dX dl dX 



(14) 



= r~'j'X + r~^j'X 



(15) 



The elasticity of the wires must be taken into account in order to increase motion control 
accuracy. The longitudinal deformation of a wire can be given by 



Then the stiffness of a wire is L 



EA 



' EA 
EA 



(16) 



, To summarize. Equations (11), (12), (13), 



/o, /,(1-A/,//,) 
(14) and (15) represent the dynamic model of the suspension system in wind tunnel. 



5. Control scheme 

The dynamic model of the flexible suspension system in wind tunnels is a highly-coupled 
and nonlinear system, and the actuation redundancy makes the system over-restrained. In 
designing the control scheme, it is necessary to decouple and linearize the dynamic model. 
A computed torque controller in joint space is employed for the flexible suspension system 
in wind tunnels. Because the actuation redundancy introduces multiple wire tension 
solutions, an optimal tension distribution method is applied to obtain certain acceptable 
solutions. When the air flow passes through the aircraft model in wind tunnel tests, a wind 
pressure will be exerted on the aircraft model. According to the aerodynamic theory, drag 
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force in along-wind direction, lift force in crosswind direction and pitching moment are 
applied on the model under the condition that the wind load is symmetrical. The wind load 
make the model fail to keep the desired position and orientation. Thus, it is challenging to 
obtain the accurate mapping relation between the measured value of the aerodynamics and 
the position and orientation of the craft model. Consequently, it is desired to enhance the 
stiffness in the three directions by commanding the wire tensions. Further more, by 
introducing the dynamical compensation on the basis of conventional PD control, the 
revised feedforward PD control law based on the stiffness enhancement principles. As 
shown in Fig. 3, the control law consisting of inverse dynamics feedforward and feedback 
loop is employed to control the driving torque of the actuators. 






^,,^,X 



Inverse 
kinematics 



Limits of dnving 

torque and 

tension 



— Ih^h 



K, 






"V Tension 
optimization 



Driving 
system 



0,0. 



Fig. 3. Control scheme for wire-driven parallel support system for wind tunnels 
The revised PD feedforward control law is 

T = A0,+C0,+rT,+K(0, -0) + K,(0, -0) 



(17) 



where, Kp , K^ are feedback gain matrices. Td is the desired tension. Error e = 0-0^. If the 
desired angular velocity 0^ , angular acceleration 0^ and tension T^ are all boundary 
values, Eq.(17) can make e and e exponentially converge to the closed sphere of radius n. 
Provided the desired trajectory Xd of the aircraft model, the desired angle, angular velocity 
and angular acceleration of the driving motors can be solved for by using inverse kinematics 
and the elastic deformation Eq.(16). 



T^=T^+Null(J)A 



(18) 



T, =r{M{X,)X, +N{X,X)^d -K) 



(19) 



where T^ is the minimal norm solution. Null is the null space vector. The restrictions for 
sinele wire and the torque of the motors are t >t>t-, where C^v is the maximum 

i-J i- maX I mm llid-A 

permissive tension of the wire, and t^^^ is the minimum tension of the wire in case of the 
pseudo drag, t^^^ and r^^^ are the maximum and minimum output torque, respectively. 
Further, we can obtain 
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(p>T>T] (20) 

where, ^. =min{^^,,(r^, -«.6;.^ -c,-4j)/^} , li = max{^^^,(r^^^ -a-^^ -^,-4^)/^} , 

For redundant driving system, an optimization is needed to solve for the tensions of the 
wires. Generally, the tension minimization principle is used in the optimization (Mtiller, 
2005). That is to say t/WT^ and W are the summing weights factor matrices. However, 
since this parallel robot is applied in wind tunnel, the aircraft model tends to deviate from 
desired position and orientation and results in experimental errors. Consequently, in order 
to obtain more precise experimental data, it is necessary to enhance the stiffness by 
adjusting wire tension. Taking account of constrains of the motor output and wire material 
properties, the objective of the optimization is to maximize the stiffness weight-sum in the 
three principal directions of forces or torques. 

find A = [\ A2] 

max sum(K^/^ ,K,/^ ,K^/^) (21) 

S.L (p(0,,O,)-f,>Null(J)A>7j(Oj,O,)-f, 

Given the desired trajectory of the aircraft model X^, T^ , (p(6^ , 9^ ) and r]{0^ , 0^ ) can be 
solved from Eq. (19) and (20). The translational stiffness in X direction is K^=K^^ = 

8 ^2 8 1 _ 

kYu-^ + Z — ^^d,i + \Null{J\^ + A^NulliJ^^) ■ ^^^* i^ ^ direction is 

8 ^2 8 1 _ 

^zz = ^3,3 = ^'Z^^ + Z — (^d,i + \Null{J).^^ + A2Null(J\2) ■ The stiffness in the pitching 

^([%X][(X''''^)X]-[(X'''^)X]['^X])2 2 - 

direction is ^^^ = ^55 = ^ vi__^_iu_^ ^) \ v\ o ^Lil_^JZ^(j^ . +X,Null{J\, 

i=\ h 

+A2Null(J)^^2) +^ Z ^ ^^ • The subscript indicates the row element of a vector or the 

element of a matrix. For the dimensionally generalized Kj^ , K^z and Knn , the objective 

function is derived by weighting sum. And the weight sum factors gi, gi and gs are 
determined according to the desired trajectory and index of the experiment. The 
optimization objective is the linear function of X^ and X2 , while the constraint function 

constitutes two-dimensional convex set of Xy and X^ . 

For this kind of linear program problem, the simplex search method is generally employed 

to solve the solution. But the solving course is very time consuming. To improve the 

computation speed, a new algorithm is designed as follow. 

Step 1: Determining the initial solution. There are sixteen linear inequality constrains in 

Eq.(21). Any three can be picked out and converted into equality constrains. Then the three 

line equation related to \ and X2 from the geometry point of view is obtained. The three 

intersection points of the three lines can be solved. Then whether the three intersection points 
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satisfy the rest thirteen inequalities is checked. If so, the three intersection points generate the 
initial solution by forming a convex combination. If not, another selection is needed. 
Step 2: Determining the searching direction. Taking the initial solution obtained in Step 1 as 
the start point. Along the gradient and negative gradient direction of the objective function 
forward search step are conducted, respectively. Then the new two candidates are evaluated 
with respect to the objective function. The direction relating to the better candidate is taken 
as the searching direction. So, this optimization becomes a one dimensional optimization. 
Step 3: Along the searching direction search is conducted forward with larger step until 
exceeds the feasible region. Then the dichotomy is used between the outer and inner points 
of the feasible region until the optimal point on the boundary of the feasible region is 
obtained. 

Step 4: In order to maintain the continuity of the wire tension, a judgment of the tension 
vector is conducted. In which, whether the tension Tj of the current position and orientation 
and Tf_i of the previous position and orientation satisfy \\T- - T-_^ ||^ ^ ^ is judged, where s is 
the threshold. If it is satisfied, the optimal solution is obtained. If it is not satisfied, starting 
with the current solution, along with the positive searching direction the optimization is 
moved back to the feasible region and the inferior solution is obtained. And the 
optimization goes to Step 1. 

6. Simulated results 

In order to validate the proposed algorithm in this research, simulations aiming at the 
revised PD feedforward controller based on the stiffness enhancement are conducted. 
Moreover, a comparison between that of a revised PD controller based on tension 
minimization is carried out. 

The position of the joints and pulley of the robot is shown in Table 1. The wire is chosen 
from reference (Zheng, 2004), which is made of extra strong polyethylene fibre. The 
diameter A = 1mm and the Young's modulus is E = 120GPa . The unit stiffness of the wire 
is ^=94247N. The maximum elastic tension is t^^^=l200^. The preset minimum 
pretension is ^j^-j^=10N. The rating output torque of the motors is r^^^^ =15.8Nm, 
Tjj^-j^ = -15.8N-m. The equivalent moment of inertia on the shaft of the motors is 
7.52xl0~'^kgm^ . The equivalent viscosity coefficient on the shaft of the motors is 
LSSxlO'^^N-m-s . The radius of the wrench is r^. = 0.04m . The ratio of the reducer is 4:1. 
The scale model is the 1/18 wooden aircraft model referred to in (Liu et al, 2005). The 
aircraft has a length of 713mm and wing width of 510mm. The height is 107mm and the 
weight is 10.5N. In the local frame, the inertial tensor is 



L 



In the experiment, the stable wind with the velocity of 30m/ s is applied. Considering the 
real-time measured data, the equivalent of load force of the wind is generated in MATLAB. 
The position of the aircraft is ^P^, = (O 420) mm and its pitch angle varies according to the 
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following parameter. The desired angle trajectory is p^= 71 1 2>^-copt and angular velocity is 
&>^ = ;r/45 rad/s , < f < 4s . In order to achieve the precision of the positioning and pitching, 
the weight factor is g^ = 0.3 , g2 = 0.3 and g3 = 0.4 . The simulation is conducted employing 
the fourth order Runge-Kutta method. The sampling time is 0.01s. The PD parameters are 
determined by both extension critical proportion and manual adjusting method. 
^^ =J/ag(l 0,10,5.5,5.5,5.5,5.5,5.5,5.5) , Kd = diag (0.25,0.25,0.15,0.15, 0.15,0.15,0.15,0.15). 



Indication 


Coordinates (in mm) 


^'Pli^'PSr^'P?) 


(-438 0)T 


o'Pli^'P^r^'Ps) 


(275 0)T 


O'p^ 


(0 -255 0)T 


^'Pe 


(0 255 0)T 


o'C 


(-25 0)T 


oBi{oB2) 


(0 0)T 


oB3{oB,) 


(0 -605 420)T 


'Bs{oBe) 


(0 840)T 


'Bvi^Bs) 


(0 605 420)T 



Table 1. Location of the joints and pulleys 

Figs. 4 and 5 show the wire tension based on the stiffness optimization principle and the 
wind load variation in three directions in the process of adjusting the orientation of the 
model. Figs 6 and 7 show the wire tension based on the minimum tension principle and the 
variation of the stiffness in three directions. Fig 8 shows the actual variation curve of the 
pitch angle in the cases of the two principles. Fig 9 makes a comparison of the positioning 
error of the aircraft model in X and Z directions. 

As the figures show, the wire tension with the stiffness optimization principle varies evenly. 
The pitching stiffness Knn ranges from 8300 to 8700 Nm/rad. The pitch angle error is less 
than 0.039 rad, and its RMS is 0.0157 rad. The positioning error in X direction is less than 
0.0111m, and its RMS is 0.0057m. The positioning error is less than 0.013m, and its RMS is 
0.0041m. 

Under the condition of minimum wire tension principle, the wire tension is small and varies 
evenly. The pitching stiffness ranges from 7600 to 8000 Nm/rad. The pitching angle error is 
less than 0.0619 rad, and its RMS is 0.0229rad. The positioning error in X direction is less 
than 0.0185m, and its RMS is 0.0076m. The positioning error in Z direction is less than 
0.0178m and the RMS is 0.005m. 

Though the wire tension based on optimal stiffness principle tends to be large, compared 
with that of minimum tension principle the pitching stiffness increases from about 7600- 
8000 Nm/rad to about 8300 - 8700Nm/rad. The RMS of the pitch angle error decreases by 
31.44%. The RMS of the positioning error in X direction decreases by 25%, and that in Z 
direction decreases by 18. The control precision has been improved drastically. 
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Fig. 4. Optimum tension distribution based on the stiffness enhancement criteria 




Fig. 5. Stiffness values obtained by the stiffness enhancement criteria 
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Fig. 6. Optimum tension distribution based on the minimum tension criteria 




Fig. 7. Stiffness values obtained by the minimum tension criteria 
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Fig. 8. Pitch angle vs time 
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7. Conclusions 

Firstly, the stiffness of the six-degree-of-freedom redundant wire driven parallel 
manipulator is dealt with in this paper. The analytical expression of the stiffness is 
developed, in which the stiffness consists of two parts. The former part is related to the wire 
tension, while the latter one depends mainly on both the geometry distribution of the wires 
and the orientation of the end-effector. 

Secondly, the dynamical models of the aircraft and the driving system are deduced, 
respectively. Considering the motor output and wire material properties, the wire tension 
optimization is conducted in order to improve the stiffness in three principal directions. This 
method solves the indefinite problem of the wires tension introduced by the redundancy. 
Thirdly, aiming at the nonlinearity, strong coupling and air current loaded environment of 
the wire driven system, the revised PD feed forward controller in joint space based on 
stiffness enhancement principle is developed. Compared with the revised PD feed forward 
controller based on minimum wire tension principle, the control scheme proposed in this 
paper improves the dynamical positioning precision of the aircraft. 
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1. Introduction 

In September 2007, a Plasma Wind Tunnel (PWT) Test was performed in the CIRA 
SCIROCCO facility on the FLPP Snecma Propulsion Solide (SPS) Thermal Protection System 
(TPS) demonstrator (Barreteau et aL, 2008). Aim of the test was to verify, in a space 
qualifying environment, the behaviour of a large assembly constituted by Ceramic Matrix 
Composite (CMC) shingles, one curved and two flat panels, the same elements which will 
be part of the next ESA Intermediate Experimental Vehicle (IX V) thermal protection system. 
The focus of this chapter is the description of the CFD activities carried out in order to 
realize and support the plasma wind tunnel test, both in the phase of test definition and for 
the post test analysis. 

During the pre-test CFD activity the test condition, previously defined by a simplified two 
dimensional methodology (Rufolo et aL, 2008), has been verified by means of three 
dimensional simulations, and the final PWT test condition has been defined. Then, the post- 
test CFD rebuilding activity has allowed the analysis of results and the comparison with 
experimental measurements. 

In addition, an assessment of the uncertainty level related to the satisfaction of the test 
requirements, in terms of heat flux and pressure to be realized over the test article, has been 
provided by analyzing the sources of error linked to both design and testing phases. 

2. Test requirements 

The test article is an assembly of CMC TPS elements: two flat panels located at 45 degrees 

with respect to the plasma flow and a curved panel which constitutes the model leading 

edge. 

The test article configuration and its dimensions are represented in Fig. 1. Each portion of 

TPS to be tested (in white in figure) is separated by the other ones by a gap (1.5 mm in depth 

and 3 mm in width), in such a way to form a ''T-gap'' configuration. 

The initial test design phase had been carried out in order to answer to the following 

customer requirements: 

• cold-wall (Tw=300 K), fully catalytic heat flux of 320 kW/m2 ± 10% at the beginning of 
the flat panels; 

• constant wall pressure of 25 mbara maintained during the test on the two flat panels 
surface. 
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Moreover, the test article leading edge (curved panel) should not have been submitted, 

during the test, to a heat flux exceeding the value of 700 kW/m^. 

At the end of the test design activity (Rufolo et aL, 2008), the PWT operating condition 

allowing the complete fulfilment of above requirements had resulted to be, in terms of 

facility reservoir conditions, Po=5.2 bara and Ho=16.7 MJ/kg, realized with the SCIROCCO 

conical nozzle D, characterised by a 1.15 m nozzle exit diameter, and with the model 

stagnation point located 0.35 m downstream of the nozzle exit section. 

The achievement of the desired operating condition (Po, Ho) in test chamber is assured by 

the measurements of stagnation heat flux and pressure on a water cooled copper probe. The 

stagnation values corresponding to the reservoir conditions above, and determined by CFD 

during the test design phase, were: Ps=36.15 mbara and Qs=2070 kW/m^. 

A complete description of the SCIROCCO facility is given in the following section. 





Fig. 1. Test article geometry 



3. SCIROCCO plasma wind tunnel 

The CIRA SCIROCCO Plasma Wind Tunnel (Marini et al, 2002 and De Filippis et al, 2003) 
is devoted to aerothermodynamic tests on components of aerospace vehicles; its primary 
mission is to simulate the thermo-fluid-dynamic conditions suffered by full scale Thermal 
Protection System (TPS) of space vehicles re-entering the Earth atmosphere. 
SCIROCCO is a large size facility (see Fig.2), whose hypersonic jet impacts the test article 
with a diameter size up to 2 m and reaches Mach number values up to 11. The jet is then 
collected by a long diffuser (50 m) and cooled by an heat exchanger. Seventy MW electrical 
power is used to heat the compressed air that expands along a converging-diverging conical 
nozzle. Four different nozzle exit diameters are available: 0.9, 1.15, 1.35 and 1.95 m, 
respectively named C, D, E and F. 

The overall performance of SCIROCCO in terms of reservoir conditions is the following: 
total pressure (Po) varies from 1 to 17 bar and total enthalpy (Ho) varies from 2.5 to 45 
MJ/kg. Enthalpy values between 2.5 and 10 MJ/kg are obtained using a plenum chamber 
between the arc heater column exit and the nozzle inlet converging part, which allows 
transverse injection of high pressure ambient air to reduce the flow total enthalpy. 
The energetic heart of the facility is the segmented constricted arc heater, a column with a 
maximum length of 5.5 m and a bore diameter of 0.11 m. At the extremities of this column 
there are the cathode and the anode between which the electrical arc is generated. A power 
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Fig. 2. SCIROCCO Plasma Wind Tunnel aerial view 

supply feeds the electrical DC power to the electrodes for the discharge. A compressed air 
supply distributes dry compressed air to the various segments of the arc heater column, 
being able to supply a mass flow rate ranging from 0.1 to 3.5 kg/s, heated up to 10000 K. 
The last important subsystem of SCIROCCO is the vacuum system, which generates the 
vacuum conditions in test chamber required by each test. The system consists of ejectors that 
make use of high pressure water steam as motor fluid (28.5 barg and 250 °C). 
Facility theoretical performance map in terms of reservoir conditions produced by the arc 
heater is shown in Fig. 3. 
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Fig. 3. Arc heater theoretical performance map 

The achievement of the desired operating conditions (Po, Ho) in test chamber before the 
insertion of the model is assured by the measurements of stagnation pressure (Ps) and 
stagnation heat flux (Qs) radial profiles on a lOOmm-diameter hemi-spherical calibration 
probe, made of copper and water cooled, at a section 0.375 m downstream of the conical 
nozzle exit section, by means of high precision pressure transducers and Gardon-Gage heat 
flux sensors, respectively. Facility regulations (mass flow, current) are tuned in order to 
measure on the calibration probe a certain couple of values (Ps, Qs) which corresponds to the 
desired set point in terms of the couple (Po, Ho). 



60 Wind Tunnels 

4. Numerical methodology 

4.1 Numerical tool 

All the three-dimensional numerical computations presented in this chapter have been 
performed by using the CIRA CFD code H3NS. 

H3NS is a structured multi-block finite volume solver that allows for the treatment of a wide 
range of compressible fluid dynamic problems, and has been widely validated in the past 
(Ranuzzi & Borreca, 2006), (Di Clemente, 2008). 

It solves the full Navier-Stokes equations for a real gas in thermal and chemical non- 
equilibrium conditions. The governing equations, written in conservation form, are 
discretized by using a finite volume technique with a centred formulation; the inviscid 
fluxes are computed by means of a Flux Difference Splitting (FDS) Riemann solver, with a 
second order ENO reconstruction of interface values, whereas viscous fluxes are calculated 
by central differencing, i.e. computing the gradients of flow variables at cell interfaces by 
means of the Gauss theorem. Time integration is performed with an explicit Euler forward 
algorithm and a Local Time Stepping formulation, coupled with a point-implicit evaluation 
of chemistry and vibrational source terms. 

In the case of thermo-chemical non equilibrium flows the fluid is treated as a mixture of 
perfect gases. The chemical model for air is due to Park (Park, 1989) and it is characterized by 
17 reactions between the five species (O, N, NO, O2, N2), neglecting the presence of inert gas or 
water in the air. The energy exchange between vibrational and translational modes is modeled 
with the classical Landau-Teller non-equilibrium equation, with relaxation times taken from 
the Millikan- White theory (Millikan & White, 1963) modified by Park (Park & Lee, 1993). The 
viscosity of the single species is evaluated by a fit of collision integrals calculated by Yun and 
Mason (Yun & Mason, 1962); the thermal conductivity is calculated by means of the Eucken 
law; the viscosity and thermal conductivity of the gas mixture are then calculated with the 
semi-empirical Wilke's formulae. The diffusion of the multi-component gas is computed 
through a sum rule of the binary diffusivities of each couple of species (Kee et al., 1983). 
Transport coefficients, assuming ideal gas, are derived from Sutherland's law. 
Several models for the treatment of finite rate catalysis are implemented both considering a 
constant recombination coefficient and ad hoc developed model for TPS materials (e.g. Di 
Benedetto & Bruno, 2010). 

4.2 Three-dimensional computational grid 

The three-dimensional computational grid around the test article has been generated by 

means of the commercial software ANSYS ICEMCFD®. 

Grid, composed of hexahedral elements, has been generated for half model using a 

multiblock approach, and has been stretched normally to wall surfaces in order to properly 

predict the different boundary layers developing around the geometric configuration. The 

topology of the grid has been created in order to accurately define all the geometric details 

of test article and obtained by using a certain number of O-grids (Fig. 4) for the block 

decomposition. 

The computational grid on the full test article is shown in Fig. 5, while an enlargement of 

top frame is depicted in Fig. 6. 

Within the main O-grid containing the body, two O-grids have been generated around the 

curved and flat panel, respectively; in this way it is possible to keep down the overall 

number of grid points still preserving a good discretization of the gaps (see Fig. 7). 
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It must be underlined that these gaps have been modelled with sharp edges (a measure of 

local curvature radii was not available), therefore results in terms of heat flux peaks are 

conservative. 

Moreover, the bow shock wave surface has been properly fitted. In order to minimize the 

numerical instabilities that propagate from the shock wave towards the stagnation region 

(the ''carbuncle'' phenomenon), it is important to align as much as possible the grid lines to 

the shock. 

Grid characteristics are listed in Tab. 1, Afimin being the minimum spacing normal to the wall 

at the stagnation point and AR the corresponding aspect ratio. Three grid levels have been 

adopted, in order to assure grid convergence of results, as it will be shown in Section 5.2. 



Grid Level 


Cells 


^«™« 


AR 


coarse 


32468 


7'10-V 


2500 


medium 


259744 


2'10"^m 


4050 


fine 


2077952 


MO"V 


4410 



Table 1. Computational grid characteristics 




Fig. 4. Block decomposition 




Fig. 5. Three dimensional computational grid 
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Fig. 6. Detail of the top frame 




Fig. 7. Gap between panels and frame (left) and T gap (right) 

5. Pre Test CFD activity 

In this section CFD three dimensional results of the flow field computed around the test 
article are shown and deeply analyzed for the PWT condition selected during the test design 
phase (Rufolo et al, 2008), i.e. Ps=36.15 mbara and Qs=2070 kW/m2. Subsequently, grid 
convergence of results will be shown in Section 5.2, and an assessment of the uncertainty 
level linked to both design and testing phases will be presented in Section 5.3. 



5.1 Three-dimensional results and test requirements verification 

Three-dimensional computations on the full test article configuration have been performed 

with the aim at verifying the test requirements fulfilment with the PWT condition defined. 

Moreover, information about flow features (presence of vortex structures, separation and 

reattachment lines, overheatings induced by the gaps, etc.) and spanwise effects will be 

given in the following, in order to exactly account for the overheatings predicted on the 

lateral parts of the CMC panels. 

The computation has been performed for half model and in the hypothesis of cold (Tw=300 

K) and fully catalytic wall, as requested by SPS at the end of the test design phase. 

Mach number and pressure contour maps are shown in Fig. 8. The shape of the bow shock 

around the model is clearly predicted as well as the stagnation pressure region (on the 

curved panel), the constant pressure region on the model flat panel and the strong 
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expansions occurring in correspondence of the roundings, either on the top frame either on 
the lateral fairings. 
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Fig. 8. Mach number around the model (left) and pressure contour map (right) 

The first verification has concerned with the possibility of wind tunnel blockage occurrence 

due to the large size and bluntness of the FLPP-SPS model. As shown in Fig. 9, where the 

computed two-dimensional and three-dimensional bow shock shapes in the model centre 

plane are reported, evident finite span effects are present for this test article which make the 

bow shock closer to the TPS demonstrator with respect to the design solution. 

The reason is the spanwise flow induced by the strong transversal pressure gradient, due to 

the 45 deg inclination of the panels with respect to the free stream. 

Fig. 10 shows the model with its bow shock wave inside the test chamber and in front of the 

diffuser entrance, at the position of 0.35 m downstream of the nozzle exit section. It is 

evident that the bow shock wave is fully swallowed by the diffuser pick-up. 

This occurrence constitutes a necessary condition to be verified in order to exclude the risk 

of wind tunnel blockage. 




L. 



Fig. 9. Bow shock in the symmetry plane 
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Fig. 10. Side (left) and front (right) view of the model with its bow shock ahead the diffuser 
entrance 

Fig. 11 shows the heat flux distribution predicted on the full model together with the skin- 
friction lines pattern (the solution on half model has been mirrored with respect to the 
symmetry plane). 

The stagnation line on the curved panel and the local maximum values of heat flux (less 
than 1 MW/m2) at the roundings of the lateral fairings of the curved panel can be clearly 
observed in the same figure, as well as the strong three-dimensionality of the flow over this 
model, that also affects the region close to the symmetry plane, where test requirements 
have been defined and matched in the test design activity (Rufolo et al., 2008). 
An enlargement of the model top frame is reported in Fig. 12, where the skin friction lines 
are coloured depending on the local shear stress value. The local maxima of shear stress are 
predicted at the shoulder of the top frame and at the roundings of the lateral fairings, as 
expected, due to the turning of the flow with associated boundary layer thinning. 
A large separated area (with negative values of shear stress) is clearly visible on the top 
frame caused by the local shock wave boundary layer interaction, with a nearly straight 
separation line and a highly distorted attachment line; the extent of the separated flow area 
increases at the extremities due to the inlet of the flow turning around the model. 
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Fig. 11. Heat flux contour map with skin-friction lines 
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Fig. 12. Enlargement of the model top frame; skin-friction lines coloured by the shear stress 

The lower frame heat flux contour map and the related skin friction lines are reported in Fig. 
13, showing a nearly two-dimensional recirculation induced by the presence of the step, 
with maximum heat flux values ranging from 45 kW/m^ in the central lower frame area to 
90 kW/m2 at the edges, where flow recirculation disappears due to the particular 
transversal shape of the model bottom part. 

The flow inside the longitudinal gap existing between the two flat panels, and inside the 
transversal gap between the full span curved panel and the two flat panels (T-gap 
structure), is described in detail from Fig. 14 to Fig. 16. A flow recirculation is predicted 
inside the longitudinal gap (see Fig. 14), with a complex vortex pattern in the ''T-gap'' region 
(see Fig. 15). The vortex flow inside the transversal gap is characterized by a strong 
spanwise velocity component, that increases moving towards the edge, a inner vortex at the 
base of the panel and an attachment line at the front edge of the panel, where very high heat 
flux values (~1 MW/m^) are predicted in a very small region. 

Fig. 16 describes the exit of the transversal gap flow into the external flow developing on the 
lateral fairing. The interaction of the two streams causes a rapid turning of the transversal 
gap flow with the formation of a local saddle point. It should be also underlined the 
presence of a inner vortex developing parallel to the junction between the flat panel and the 
lateral fairing, and the presence of an attachment line (the same already seen in Fig. 15) at 
the front edge of the flat panel, which corresponds to a region of high heat flux, with a 
maximum in the top corner of about 1.6 MW/m^ but localized in a very small region (0.0002 
m depth). 

In order to verify test requirements in terms of heat flux and pressure at the beginning of the 
flat panel, and to properly evaluate spanwise and viscous effects, the longitudinal and 
transversal distributions along the slices indicated in Fig. 17 have been analyzed. 
Results in terms of heat flux are reported in Fig. 18 and Fig. 19, showing transversal and 
longitudinal distributions, respectively, these latter ones compared to the two-dimensional 
results of test design activity (Rufolo et al., 2008). 
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Fig. 13. Heat flux contour map with skin-friction lines; model bottom frame 




Fig. 14. Re-circulating region; longitudinal gap 
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Fig. 15. T-gap; heat flux contour map with skin-friction lines 
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Fig. 16. Exit of transversal gap flow. Heat flux contour map and skin-friction lines 




Fig. 17. Longitudinal and transversal slices 

The increase of heat flux predicted on the flat panel is due either to spanwise effects either to 
the presence of gaps (longitudinal and transversal) and steps (lateral side), as clearly shown 
in Fig. 18. At the flat panel leading edge three-dimensional CFD simulation yields a 28% 
increase (450 kW/m^) of predicted heat flux, both 5mm from the centreplane (Z=0.005m) 
and 5mm from the lateral edge (Z=0.195m), and it is nearly 350 kW/m^ in-between. 
Downstream along the panel the predicted heat flux is closer to the test requirement, while 
localized high heat flux peaks are present in correspondence of gaps and steps. 
Transversal and longitudinal wall pressure distributions are shown in Fig. 20 and Fig. 21, 
respectively. Pressure is not affected by spanwise effects from the qualitative point of view 
(the transversal distributions remain two-dimensional for most of the half panel span), but a 
quantitative reduction of 17% of maximum pressure on the flat panel centreplane is 
predicted (2070 Pa instead of 2500 Pa). 
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Fig. 18. Transversal heat flux distributions 
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Fig. 19. Longitudinal heat flux distributions; comparison with 2D distribution 
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Fig. 20. Transversal wall pressure distributions 
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Fig. 21. Longitudinal wall pressure distributions; comparison with 2D distribution 

5.2 Grid convergence of results 

Grid convergence study is the most common and reliable technique for the quantification of 

numerical uncertainty (Roache, 1998) related to spatial discretization. It has been carried out 

for the three-dimensional pre-test computation by using the different grid levels indicated in 

Tab. 1. 

Temporal convergence of the solutions has been obtained on all the grid levels. 

Grid convergence of results has been evaluated in correspondence of the same points used 

in the design phase for monitoring the test requirements matching, i.e. the beginning of flat 

panel for the heat flux and the point of maximum value for the pressure on the flat panel, 

both taken at the centreline. In the three-dimensional case, these control points have been 

selected in the spanwise direction in order to be close to the symmetry plane, but sufficiently 

far from the region affected by the presence of the longitudinal gap; their coordinates are 

reported in Tab. 2. 

Q* and P* indicate the values of heat flux and pressure in the selected points. 



z=0.07 m 


X 

(for Q evaluation) 
-0.172 m 


X 

(for P evaluation) 
-0.156 m 



Table 2. Coordinates of the points selected for the grid convergence study 



GRID 


N 


N-1/3 


Q*(W/m^) 


P*(Pa) 


coarse 


32468 


0.0313 


132675.69 


1959.30 


medium 


259744 


0.0157 


335118.84 


2024.86 


fine 


2077952 


0.0078 


349148.53 


2044.60 


Rich.Extrap. 


inf. 





353825.09 


2051.19 



Table 3. Q* and P* values at the selected points for the three grid levels and Richardson 
Extrapolation 
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The computed Q* and P* values are reported in Tab. 3 for the three grid levels, together with 
the Richardson Extrapolation value. This latter is an estimation of the ''continuum value'' 
(i.e., the value at zero grid spacing), obtained from a series of discrete values, and it is 
defined in the following way: 

//.=o^/i+4^ (1) 

r^ -1 

where: fh^o is the value at zero grid spacing; /i and/2 are the values computed on two grids, /i 
being the finer one; p is the order of the solution {p=2 for this case); r is the grid refinement 
ratio: 



(2) 

Ni and N2 being the numbers of cells of the grids 1 and 2, respectively. In the following, N 
will be used to indicate the total number of cells of a grid level, while (l/N)-'^/^ is a parameter 
that represents adequately the grid resolution. 

The difference between the values /i and fh=0 is one of the error estimators. The actual 
fractional error is defined as: 

v4i=AzA^ (3) 

fh^O 

Another error estimator, the relative error, is based on the difference between /i and /2: 




/l 



(4) 



This quantity has to be corrected to take into account r and p. The estimated fractional error 
for/i is therefore defined as: 

^1=1^ (5) 

r^ -1 

Although El is based on a rational theory, it is not a bound on the error. On the contrary the 
Grid Convergence Index (GCI) provides an error band, i.e. a tolerance on the accuracy of the 
solution (Roache, 1998). The GCI on the fine grid is then defined as: 

GCI,„=J^ (6) 



fine 



("-') 



where Es is a safety factor, that is recommended to be 3.0 when comparing the results of two 
grids, and 1.25 for comparison of three grids (being this latter our case). The above defined 
error estimators have been all calculated, and are reported in Tab. 4 for Q* and P*. 
The values of heat flux (Q*) and pressure (P*) are reported in Fig. 22 for the three grid levels 
in function of the grid resolution (i.e. the parameter (1/N)-^/^) and compared with the value 
corresponding to zero grid spacing (computed by means of the Richardson extrapolation). 
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Error Indices 


Q*(W/m^) 


P*(Pa) 


eps 


0.0402 


0.0097 


E1 


0.0134 


0.0032 


GCI 


0.0167 


0.0040 


A1 


-0.0132 


-0.0032 



Table 4. Grid error indices 

These plots confirm the right trend of solution grid convergence both for heat flux and 
pressure. In fact, the difference existing between the results of the coarse grid level and the 
medium one decreases if comparing the medium level with the fine one, and the trend of 
solution is towards the Richardson extrapolated value. 

As a consequence, the Grid Convergence Index provides a level of confidence of the 
solution, therefore it can be concluded that (see Tab. 4): 

• the error committed on the heat flux value with the finer grid level should be lower 
than 1.67%; 

• the error committed on the pressure value with the finer grid level should be lower 
than 0.40%. 
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Fig. 22. Grid convergence estimation for heat flux (Q*) and pressure (P*) at the selected points 



5.3 Estimation of uncertainties 

An assessment of the uncertainty level related to test requirements fulfilment in terms of 
heat flux and pressure to be realized over the test-article is provided in this subsection, both 
for test design and test execution phases (Rufolo et al., 2008). The high complexity of 
involved phenomena together with the heterogeneous character of the different error 
sources make it impossible to give a rigorous definition and quantification of the error, but 
only a simplified estimation can be pursued. 

Fig. 23 reports the entire process of numerical test design and test execution: during the design 
phase, starting from test requirements, a CFD aided activity is carried out in order to derive 
the proper settings for the heat flux (Qs) and pressure (Ps) over the PWT calibration probe; in 
the testing phase the facility driving parameters (mass flow and arc current) are tuned in order 
to get the desired couple (Qs, Ps) over the calibration probe, then the test is executed and with 
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Fig. 23. Numerical test design and test execution chain 

the post-test analysis it is finally possible to verify the matching of the requirements. 
Obviously, an error £f is linked to each phase of the above described chain, and all of them 
contribute in determining the difference between the original requirements and their actual 
realization. It has to be said that in the present case the requirements were expressed in terms 
of heat flux and pressure for a fully catalytic and isothermal cold wall, and this is clearly a not 
realistic hypothesis for the kind of material and type of test to be conducted. 
Moreover, during the test no heat flux direct measurements have been provided, and only 
an indirect derivation from temperature measurements can be obtained assuming radiative 
equilibrium at the wall (i.e. neglecting conduction into the material). In order to fully exploit 
measurements it is needed to associate correct values of catalytic recombination and 
emissivity coefficients, but these data have not been available during the project. 
For these reasons, being unfeasible to characterize the complete error chain, only the 
following components of the error chain will be described hereinafter (Rufolo et al., 2008): 

• how the test requirement is translated by means of CFD into PWT conditions (£2 in Fig. 23); 

• how the error in the experimental realization of the set point propagates on the 
requirements over the test-article (£3 in Fig. 23). 

The evaluation of the error £3 propagation is made by substituting the facility with its 
numerical modelling. 

The numerical setting of PWT operating conditions comes out from an iterative process in 
which the facility driving conditions (Ho, Po) are tuned in order to match the requirements in 
terms of heat flux and pressure over the model to be tested (Di Benedetto et al., 2007). The 
error related to this process is definitively negligible, in the sense that it is always possible to 
find a couple (Ho, Po) that allows to numerically satisfy the requirements whichever is the 
accuracy prescribed. At the end of this process, when the correct couple (Ho, Po) has been 
found, the simulation of the flow field around the calibration probe is carried out in order to 
find out the couple (Qs, Ps) that will be used for the test execution (Di Benedetto et al., 2007). 
The process that translates the reservoir condition (Ho, Po) in local parameters (Qs, Ps) by 
means of a numerical modelling is affected by an error, above defined as £2. 
By following the classical taxonomy adopted for CFD (AIAA, 1998) it is possible to 
recognize the following three error components for £2: 

• the Modelling Error (Chemical processes, fluid properties. Initial and Boundary 
conditions. Geometry representation. Turbulence Model); 

• the Discretization Error (Grid independence, algorithm error); 

• the Iteration Error (Convergence criterion). 
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The modelling error is by far the most complex source of uncertainty to estimate. The 
common practice (AIAA, 1998) relies on the validation of the numerical code with respect to 
experimental measurements obtained for simple test cases. Unfortunately, the experimental 
measurement it is affected by an error that, especially in the case of heat flux measurements 
for aerothermodynamic tests, can make void the validation process. 

As reported in (Ranuzzi & Borreca, 2006) a series of comparisons with existing literature 
experiments were carried out during the development and validation phase of the H3NS 
CFD code. In particular, it was decided to refer to the Hyperboloid Flare Test Case carried 
out at the F4 blow-down arc heated high enthalpy facility of the ONERA in order to find out 
an error level applicable to the present case (Rufolo et al., 2008). The freestream Mach 
number is 8.7, the total enthalpy is about 13 MJ/Kg, the wall is considered isothermal at a 
temperature of 300 K and fully catalytic. Trying to find out an estimation of the modelling 
error related to the phenomenon we are interested in (heat flux and pressure along the test- 
article flat panel), it is possible to extract the average percentage error for the measurement 
stations located in the mid part of the hyperboloid and ahead of the flare. In this way an 
error of about 4% for heat flux and 3% for pressure is obtained. 

Another possibility for estimating the modelling error, in absence of affordable experimental 
results, is to carry out a sensitivity analysis with respect, for instance, to chemical model 
and/ or transport properties model. With respect to the transport properties model, results 
obtained for the hyperboloid flare show no significant effect on pressure, while for heat flux 
the maximum deviation is about 3.1%. As for the chemical models, a dedicated analysis has 
been carried out both for the PWT calibration probe and for the SPS test-article. The four 
different chemical models implemented in H3NS (Ranuzzi & Borreca, 2006) have been 
tested: Kang-Dunn (Dunn & Kang, 1997), Park 1990 (Park, 1990), Park-Rakich (Rakich et al, 
1983) and Park 1993 (Park & Lee, 1993), this latter being the chemical model used for all the 
simulations performed in the present activity. Regarding the stagnation point of the 
calibration probe, the largest deviation occurs for the Kang-Dunn model (2.63% for heat flux 
and 0.97% for pressure). For what concerns the SPS test-article simulation, the percentage 
deviations of heat flux at the beginning of the flat panel and of maximum pressure over the 
flat panel obtained with Kang-Dunn model with respect to the Park 1993 results are 
respectively 0.38% and 3.13%. 

For what concerns the discretization error, the results of the grid convergence analysis of the 
three-dimensional simulation of the FLPP-SPS test-article, reported in Section 5.2, show that, 
with respect to an ideal zero-spacing grid, an error of 1.67% on the heat flux at the beginning 
of the flat panel and of 0.40% on the maximum pressure on the panel is committed. 
For what concerns the iteration error, it has to be said that, even if we are interested in 
achieving the steady state solution of the Navier-Stokes equations, when the flow field to be 
resolved contains features characterized by intrinsic unsteadiness (e.g. recirculation bubble, 
vortex shedding, shock wave instability), the residue of the equations does not decrease 
towards the machine precision. Despite the presence of these unsteadiness, the quantities of 
interest in our case, as the heat flux and the pressure over the flat panel, reached a steady 
state value so that the iteration error can be neglected. 

Trying to summarize. Tab. 5 reports the identified uncertainties (intended as estimation of 
the errors). The last column of the table reports the ''overall error'' obtained adding all the 
components. 

Concerning the error £3, it is needed to estimate how the experimental uncertainty on the 
measurements of heat flux and pressure over the calibration probe translates in uncertainty 
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Modeling error 






Discretization 


Iteration 


Validation 


Chem 

Model 

Sensitivity 


Transport 

Model 
Sensitivity 


Overall 
Error 


S2(Q) 


1 .7% 


~0 


4.0% 


2.6% 


3.1% 


11.5% 


e2(P) 


0.4% 


~0 


3.0% 


3.1% 


~0 


6.5% 



Table 5. Summary of identified error components 

of the requirements over the test article. This can be done only by adopting the CFD tool as 
transfer function. The error associated to the heat flux measurement of the calibration probe 
is ±90 kW/m2, while the one associated to the stagnation pressure measurement is ±1.1 
mbar (Marini et al, 2002). Starting from the values realized during the test, four couples (Qs, 
Ps) have been identified by adding and subtracting their own error both to Qs and Ps; the 
corresponding values are reported in the first two columns of Tab. 6. For each couple, the 
facility driving conditions (Ho, Po) have been found by following the iterative process 
already described (Di Benedetto et al., 2007) (columns three and four of Tab. 6), and then 
two-dimensional simulations of both probe and model (i.e. test article) have been carried out 
for each condition. The percentage errors referred to the nominal values are reported in Tab. 
6 for each of the four conditions. 





Ps 
[mbar] 


PROBE 

Qs 
[kW/m2] 


EXP 

PO 
[bar] 


HO 
[MJ/kg] 


Ps 
[mbar] 


PROB 

err% 


ECFD 

Qs 

[kW/m2] 


err% 


MODEL CFD 

POINT #4 POINT #1 
Ps err% Qs err% 
[mbar] [kW/m2] 


(Preq.Qreq) 


34.2594 


2121.82 


4.90 


17.40 


34.26 


2121.82 


23.84 


338.66 1 


[Preq,qreq+err(q)] 


34.2594 


2211.82 


4.88 


18.03 


34.25 


0.02% 


2211.56 


4.23% 


23.82 


0.08% 


353.02 


4.24% 


[Preq,qreq-err(q)] 


34.2594 


2031.82 


4.93 


16.81 


34.26 


0.01% 


2030.47 


4.31% 


23.84 


0.02% 


325.32 


3.94% 


[Preq+err(p),q,J 


35.3594 


2121.82 


5.07 


17.18 


35.35 


3.17% 


2122.05 


0.01% 


24.58 


3.12% 


339.47 


0.24% 


[Preq-err(p),qreq] 


33.1594 


2121.82 


4.73 


17.64 


33.18 


3.14% 


2125.11 


0.16% 


23.09 


3.14% 


338.39 


-0.08% 



Table 6. Influence of calibration probe measurements uncertainty on test article 
requirements 

Regarding the model, the errors were evaluated with respect to the beginning of the flat panel for 
the heat flux (Point #1 in Tab. 6) and to the maximum value over the flat panel for the pressure 
(Point #4 in Tab. 6). It can be seen that in the worst case the difference between the errors on the 
probe and on the model are limited to 0.37% for the heat flux and 0.05% for the pressure. So it can 
be stated that, within the approximation related with the numerical process, the experimental 
uncertainties on the point settings is identically transferred to test article requirements. 
In conclusion, the analysis reported above has been aimed at deriving an estimation of the 
errors £2 and £3. Obviously, the analysis cannot be considered exhaustive and especially for 
the CFD related error only a very simplified indication has been provided. As a matter of 
fact, the two errors £2, £3 can be considerer fully independent. 

At worst, for the present case the estimated overall errors are about 15% on heat flux and 
9.5% on pressure. 



6. Rebuilding CFD activity 

The FLPP-SPS TPS demonstrator plasma wind tunnel test was successfully performed on 
September 20th, 2007 simulating a 15 min re-entry trajectory in three steps characterized by 
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increasing total enthalpy level in test chamber, i.e. increasing continuously wall heat flux 

(Trifonietal.,2007). 

The test condition, which the CFD three-dimensional analysis described in the previous 

section refers to, corresponds to the second test step, defined as the ''nominaF' one. This 

latter condition has been rebuilt after the test by exploiting the calibration probe heat flux 

and pressure available measurements. 

A different hypothesis about the temperature wall condition has been made, in order to 

simulate a more realistic condition with respect to the hypothesis of cold wall of the pre-test 

CFD simulation. In particular, radiative wall temperature has been computed assuming the 

equilibrium between the convective and the radiative heat fluxes. The emissivity coefficient 

has been provided by SPS (s=0.8), while the hypothesis of fully catalytic surface has been 

maintained also in the test rebuilding CFD simulation, as also indicated by SPS. 

In Fig. 24 the CAD model (left) is compared with the model as built (right), in which there is 

no step in the bottom part. However, this difference in the test article configuration should 

involve discrepancies only on the regions closer to the bottom part of the model, therefore 

no influence is expected on the flat and curved panels. 




Fig. 24. CAD model (left) and model as built (right) 

6.1 Operating condition assessment 

The pre-test three-dimensional CFD simulation has been carried out in the PWT operating 
condition resulting from the previous CFD test design activity (Rufolo et al., 2008), whose 
results are reported in Tab. 7. 



Design Test Chamber 
Conditions 


Po (bara) 


Ho(MJAg) 


5.20 


16.70 


Calibration Probe 
Stagnation Point (CFD) 


Ps (mbara) 


Qs (kW/m2) 


36.15 


2070 



Table 7. PWT test design operating condition 

This condition has been compared, in terms of heat flux and pressure on the PWT 
hemispherical calibration probe, with that actually measured during the second step (the 
"nominaF' one) of the test. These latter values are reported in Tab. 8, together with their 
error bars (Trifoni et al., 2007). 
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In order to reproduce in the rebuilding CFD simulation the same condition realized in test 
chamber during the test in terms of total pressure and total enthalpy, the iterative procedure 
described in (Rufolo et aL, 2008) and (Di Benedetto et aL, 2007) has been applied, this time 
having as requirements the values measured on the calibration probe. 



Calibration Probe 

Stagnation Point 

(Measured) 


Ps (mbara) 


Qs (kW/m2) 


34.20±11 


2120±90 



Table 8. Values at the calibration probe stagnation point measured during the test 

Finally, the PWT operating condition obtained for the rebuilding CFD activity is 
summarized in Tab. 9. 



Rebuilding Test 
Chamber Conditions 


Po (bara) 


Ho(MJAg) 


4.90 


17.40 


Calibration Probe 
Stagnation Point (CFD) 


Ps (mbara) 


Qs (kW/m2) 


34.25 


2121 



Table 9. PWT test rebuilding operating condition 



6.2 Three-dimensional results 

The three-dimensional CFD rebuilding simulation has been performed in the PWT 
"nominaF' test condition of Tab. 9. The more realistic radiative equilibrium wall condition, 
with surface emissivity s=0.8, has been imposed instead of the cold wall. In order to 
qualitatively evaluate the actual catalysis of the CMC panels through comparison with 
temperature measurements, both fully catalytic (FC) and non catalytic (NC) wall conditions 
have been considered. 

Heat flux distribution together with the skin-friction lines pattern on the test article is shown 
in Fig. 25: heat flux on the stagnation line is about 600 kW/m^ for FC case, and it decreases 
to 200 kW/m2 for NC one. Temperature contour maps are shown in Fig. 26: in the FC case 
the local maximum values of temperature are around 2000 K on the stagnation line and 
about 2200 K on the roundings of lateral fairings of the curved panel. On the flat panels the 
predicted temperature ranges from about 1500 K (in the single panel central area) to about 
1800 K at the panel lateral edges. Temperature levels of about 1000 K are predicted on the 
lateral sides of the test article. These values are quite strongly reduced with the NC 
assumption (about 500 K on the stagnation line), due to a combined effect of the high energy 
content of the flow and the large bluntness of the test article. 

The analysis which follows refers to FC condition results only, this in order to make possible 
a comparison with the pre-test numerical findings. An enlargement of the model top frame 
with skin-friction lines coloured by shear stress value is reported in Fig. 27 (left) and 
compared with the distribution obtained in the pre-test simulation (right). The 
phenomenology and the shear stress distribution are very similar to those predicted in the 
pre-test activity, while a slightly larger separated area is observed as a consequence of the 
changed wall temperature condition. In fact, a higher surface temperature implies a 
boundary layer thickening (in particular of the subsonic region), in this way increasing the 
upstream and downstream pressure disturbance propagation. As a consequence of the 
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Fig. 25. Heat Flux contour map with skin-friction lines; FC (left), NC (right) 
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Fig. 26. Temperature contour map; FC (left), NC (right) 
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Fig. 27. Enlargement of the model top frame; skin-friction coloured by the shear stress; 
rebuilding (radiative equilibrium, left), pre-test (cold wall, right) 
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Fig. 28. T-gap heat flux contour map with skin-friction lines (left) and longitudinal gap 
recirculation (right) 
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Fig. 29. Transversal pressure (left) and heat flux (right) distributions 
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Fig. 30. Longitudinal pressure (left) and heat flux (right) distributions 
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increased temperature, an extension of the regions submitted to higher shear stress is 
observed, although the overall structure of the flow seems unchanged. 

The flow inside the T-gap is described in Fig. 28. The interaction between the transversal 
stream and the longitudinal one realizes in a saddle point and in two lateral vortices, but 
with a different flow pattern with respect to the pre-test simulation due to the effects of the 
surface temperature wall condition (see Fig. 14 and Fig. 15). The vortex flow inside the 
transversal gap is again characterized by a strong spanwise velocity component that 
increases moving towards the edge, a inner vortex at the base of the panel and an 
attachment line at the front edge of the panel. As expected, the region of high heat flux at the 
front edge of the flat panel, and in particular at the top corner, is largely reduced. 
Pressure and heat flux distributions in transversal and longitudinal directions are shown, 
respectively, in Fig. 29 and Fig. 30. The main flow features, already described in Section 5.1 
(see from Fig. 18 to Fig. 21), are all confirmed by the present test CFD rebuilding, although 
quantitative levels are different due either to the realization of a slightly different ''nominaF' 
condition, with respect to that analyzed during the pre-test CFD activity, either to the 
different surface thermal boundary condition. 

At the flat panel leading edge, CFD rebuilding simulation yields a heat flux of about 440 
kW/m2 5mm from the lateral edge (Z=0.195m), and it is slightly larger than 300 kW/m^ for 
the rest of the panel (Fig. 29-right). Downstream along the panel heat flux remains around 
300 kW/m2 apart from the lateral edge, affected by the presence of the step, where 400 
kW/m2 all along the panel are predicted (Fig. 30-right). 

Transversal and longitudinal pressure distributions over the model are reported in Fig. 29- 
left and Fig. 30-left respectively; pressure is not significantly affected by spanwise effects, 
apart from the more lateral section Z=0.195 m where a strong flow expansion occurs: 
transversal distributions remain two-dimensional for most of the half panel span, as well as 
the longitudinal ones are flat enough for 80% of the panel length. 

7. CFD/Experiments comparison 

In this section some of the experimental data collected during the FLPP-SPS demonstrator 
test in the SCIROCCO PWT (Trifoni et al., 2007) are compared to the results of the numerical 
rebuilding described in Section 6. 




Fig. 31. Test article instrumentation 
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During the test, eleven B-type thermocouples have measured the back wall temperatures of 
the CMC panels. Among these, those located on the flat panels which have correctly worked 
(F2-1, G2-1, H2-1, Hl-1, see Fig. 31) have been selected to perform comparisons with CFD 
temperature distributions. Moreover, a dual colour pyrometer (range: 1000-3000 °C) has 
been pointed to G2-1 thermocouple location and two IR thermo-cameras (s=0.8, range: 600- 
2500°C) have been used to monitor the test article during the test both from the top (flat 
panels) and from the lateral front (curved panel area). 

In Fig. 32 temperature measured by thermocouples is compared with CFD distributions 
along the two sections, indicated as slices in the figure, where thermocouples are located. 
As expected, measured temperatures lie more or less in the middle between the non catalytic 
(NC) and the fully catalytic (FC) distributions. In addition, it has to be said that the surface 
temperatures can be estimated to be about 50 °C higher than the measured back wall ones. 
In Fig. 33, the same kind of comparison is reported for the temperature measured by the 
dual colour pyrometer. A lower emissivity value of 0.68, which is a combination of the real 
emissivity value of the material and all the experimental uncertainty factors, allows to match 
pyrometer and thermal camera readings, as reported in Tab. 10 (experimental emissivity 
evaluation). Therefore, also the CFD temperatures in Fig. 33 have been properly scaled (to 
the emissivity value of £exp=0.68) in the post-processing phase, in such a way to make the 
comparison meaningful and to reproduce as much as possible the actual wall conditions. 
An attempt to derive an estimation of the CMC panels catalytic recombination coefficient 
has been done by combining the experimental results to a CFD-based correlation. Namely, 
by means of CFD two-dimensional computations with finite rate catalysis values at the wall, 
a function that relates the heat flux at a certain point of the flat panel with the recombination 



T 
pyrometer 


T 
thermocamera 


£exp 


1500 K 


1360 k 


0.68 



Table 10. Experimental emissivity evaluation 




Z(m) 

Fig. 32. Comparison between temperature CFD distributions and thermocouples 
measurements 
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Fig. 33. Comparison between temperature CFD distributions and pyrometer measurement 

coefficient y has been derived. By crossing this function with the radiative heat flux 
corresponding to the pyrometer reading, a value for y of about 0.008 has been obtained. It 
has to be remarked that this value only represents a rough estimation and it includes all the 
numerical and experimental errors. 

Finally, some qualitative comparisons of the bow shock wave shape are shown from Fig. 34 to 
Fig. 36, where the predicted flow field in the shock layer region has been overlapped to the 
images taken by the two video cameras during the test. In Fig. 34 and Fig. 35, the shock section 
extracted from CFD computation and the predicted temperature field in the shock region have 
been superimposed on a view from the top camera. The comparison shows that both shock 
shape and stand off distance predicted in the stagnation region well reproduce the actual ones. 
In Fig. 36 the predicted atomic nitrogen mass fraction is overlapped to a view from the side 
camera, showing a good agreement of predicted and actual shock shape around the entire 
model, and a significant presence of atomic nitrogen (N) around most of the curved panel. 




Fig. 34. Top view of the model during test. Comparison of predicted and actual shock shape 
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Fig. 35. Top view of the model during test. Comparison with predicted temperature 
contours 
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Fig. 36. Side view of the model during test. Comparison with predicted nitrogen 
concentration 



8. Conclusions 

This chapter has described the three-dimensional CFD activities carried out to support the 

SCIROCCO plasma wind tunnel test performed on the FLPP-SPS TPS demonstrator 

designed and manufactured by Snecma Propulsion Solide. 

After a CFD pre-test activity, during which the test point previously designed by a 

simplified two-dimensional methodology has been verified and the final PWT test condition 

frozen, the post-test phase has regarded the plasma test CFD rebuilding. 

The FLPP-SPS PWT test was performed with full success on September 20th, 2007 simulating 

a 15 minutes re-entry trajectory in three steps characterized by increasing total enthalpy 

levels in test chamber. The test condition which the present CFD three-dimensional analysis 

refers to corresponds to the second ''nominaF' step. 

This latter condition has been rebuilt by exploiting the calibration probe heat flux and 

pressure available measurements, and by applying the same iterative procedure used 
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during the test design phase, this time having as requirements the values measured on the 
calibration probe. Moreover, in order to perform more realistic simulations, radiative 
equilibrium has been imposed at the wall, whereas to qualitatively evaluate the actual CMC 
panels catalysis both FC and NC conditions have been considered. 

Similar flow features have been predicted both in the pre-test phase and the post-test 
rebuilding phase, and some meaningful comparison between CFD rebuilding results and 
experimental findings have allowed to assess the full capability of the present CFD-based 
methodology to design and properly rebuild a plasma wind tunnel test, with its own 
accuracy bounds. In addition, an approach to determine the uncertainties related to both 
design and testing phases, with respect to the satisfaction of test requirements, has been 
presented. 

Finally, a rough estimation of the catalyticity of the CMC panels under realistic re-entry 
conditions has been obtained by crossing experimental measurements and CFD results. 
An important step for future applications like the present should be to rebuild plasma wind 
tunnel tests accounting for the actual catalytic behaviour of the different parts of the test 
article. Of course, to do this the proper experimental characterization of the involved 
materials in terms of recombination coefficients as functions of temperature and pressure is 
needed. Then, once having re-tuned the CFD methodology, the approach could be directly 
applied starting from the pre-test design phase. 
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1. Introduction 

The modal decomposition of unsteady flowfields was proposed in the 1990s by several 
authors, e.g. Hall (1994) or Dowell et al. (1998). Proper Orthogonal Decomposition (POD) 
is one method that can be used in order to perform this modal decomposition; it became 
popular for aerodynamics research in the 2000s, starting with Tang et al. (2001), although it 
was first proposed for use in fluid dynamics in the 1960s by Lumley (1967). 
The basic principle of POD is the creation of a mathematical model of an unsteady flow that 
decouples the spatial from the temporal variations. A 2D flowfield described by the horizontal 
velocity u{x,y, t) and the vertical flow velocity v{x,y, t) can thus be expressed as 

M 

u{x,y,t) = u{x,y) + u'{x,y,t) = u{x,y) + ^qi{t)(puj{x,y) 

i=\ 

M 

v{x,y,t) = v{x,y)+v'{x,y,t) = v{x,y) + J^qi{t)(p^^i{x,y) (1) 

i=l 

where U{x,y) and v{x,y) are obtained by time averaging the flowfield over M time instances, 
while u'{x,y,t) and v'{x,y,t) are time-dependent fluctuations from the mean. These 
fluctuations are decomposed using M mode shapes (pu,i{X/y)/ (pvA^^y) ^^^ ^ generalized 
coordinates qi{t). For a reduced order model, the number of modes, N << M, is to be chosen 
as a compromise between model simplicity and model accuracy. The principle of the POD 
technique is to extract the most energetic modes that capture most of the unsteady flow energy. 
The POD technique has been used to decompose several types of aerodynamic flows, such as 
the flow behind a disk (Tutkun et al., 2008), the flow past a delta wing (Cipolla et al., 1998), 
the unsteady flow impinging on an aircraft tail behind a delta wing (Kim et al., 2005), the 
unsteady flow around a F-16 fighter configuration (Lie & Farhat, 2007) and others. 
It should be noted that there are two types of POD research being carried out at the moment. 
The first concerns the decomposition of flowfields observed in experiments in order to better 
understand the flow mechanisms and physics underlying these flows. The second type of 
research concerns the Reduced Order Modelling of unsteady Computational Fluid Dynamic 
(CFD) simulations or even, CFD/CSD (Computational Structural Dynamics) simulations. 
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in order to produce simplified but representative models that can be used in practical 

applications such as aircraft design. 

The work of interest here is the of first type, i.e. the experimental work. It is usually combined 

with high-speed Particle Image Velocimetry (PIV) measurements, although there are examples 

of other instrumentation being used, such as hot wire rakes (Tutkun et al., 2008). The limitation 

of all research works published on the subject is that the models around which the flowfield 

is measured are always static or rotating at constant velocity. Additionally, only one source of 

flow unsteadiness is ever considered. 

The objective of the present work is to expand the methodology of the application of POD to 

experimental flowfields. There are two aspects to this expansion: 

1. Allow the models to oscillate. The source of the unsteadiness will then be the movement 
of the model, as well as any unsteadiness due to flow separation. 

2. Study the interaction between the different sources of unsteadiness. In particular observe 
how the modes generated by one source of unsteadiness interact with the modes generated 
by the other. Determine if it is possible to separate the structural from the aerodynamic 
sources of unsteadiness. 

2. Basics of Proper Orthogonal flow decomposition 

Observation of an unsteady flow by PIV will, in general, yield M shapshots of a 2D section 
of the flowfield at times f^, . . . , f^- These snapshots will usually feature information on the 
u{x,y,t) and v{x,y,t) velocity vectors although other information can also be obtained (e.g. 
vorticity). The velocity vectors will be available on a spatial grid of size Uy x nj, i.e. there 
will be fly gridpoints in the y direction with spacing Sy and Ux in the x direction with spacing 
Sx. Therefore, u{x,y,t) and v{x,y,t), will be available in discrete form, i.e. in the form of 
Hy X Hx X M real arrays. 
The time-averaged flow is represented by (u{x,y),v{x,y)), where 

■y M 

u{x,y) = j^Y^u{x,yJi) 

^ M 

and the unsteady velocity components are obtained simply from 

u'{x,y,t) = u{x,y,t) — u{x,y) 

v'{x,y,t) = v{x,y,t)-v{x,y) (2) 

The Proper Orthogonal Decomposition procedure is then applied on the data matrix C, the 
auto-correlation matrix of the total energy in the flow at every instance in time. For a 
continuous flow, 

C{h,t2) = —jj{u'{x,y,ti)u'{x,y,t2) + v'{x,y,h)v'{x,y,t2))dxdy (3) 

For a discrete flow, the integrals become summations. Using trapezoidal integration. 
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1 "y-l 
^ii =mL (Gij.k + G.-,y,fc+i)^y/2 (4) 

^^^ k=l 

where C^ y is the element in the zth line and ;th column of C, G^y^ = Efi7 (^ijXl + 

^z,7',fc,/+l)^^/2/ ^i,i,k,l ~ \^'kli^'kli~^^'kli^'kli] ^^^ *^^ notation u'^n is shorthand for 

u' {x]^,yi, ti). Higher order integration schemes can also be used. 

The Proper Orthogonal Decomposition process requires the solution of the eigenvalue 

problem 

CA = AA (5) 

where A are the matrix of eigenvectors of the C matrix and A are its eigenvalues. If the 
eigenvectors are normalized in the form a^ / v^A^M, where a^- is the zth column of A, then they 
will form an orthonormal basis. The mode shapes ^^i and ^^i can be constructed from 



1 



M 



^ M 

(pvA^^y) = -TrT^ E ^'{^^y^tm)am,i (6) 

where a^ ^ is the mth element of the z'th eigenvector of C. The mode shapes are only 
functions of space but can be used to describe the unsteady flowfield when combined with 
the generalized coordinates qi{t), which can be obtained from 

^z(0 = j j {u\x,yj)(puj{x,y)^v\x,yj)(p^j{x,y)) dxdy (7) 

or from the discrete version of this equation. 

There are M eigenvalues and hence M sets of mode shapes and generalized coordinates. 

However, the aim of POD is to create a reduced order model, using only the first N modes 

that contain most of the fluctuating flow energy. To this end, the quantity A^/ Y^f^i A/ can 

be inspected, assuming that A/ is ordered from highest to lowest eigenvalue. Kim et al. (2005) 

note that it should be Ylfti ^i = 1/ however this depends on the scaling used by the eigenvalue 

estimator, so looking at the ratio A// T.f^i A/ is safer. If the first N eigenvalues are chosen, that 

have ratios higher than, say, 10% then a N mode model will be created. 

Finally, the N-mode approximation of the complete velocity field can be reconstructed from 

the N retained modes using equation 1 

N 
w*(x,y,f) = u{x,y) + J^qi{t)(l)uj{x,y) 

i=l 

N 

p*(x,y,0 = v{x,y)^J^qi{t)(l)j,j{x,y) 

i=l 

The POD procedure described above will be applied (with some modifications) to 
experimentally observed unsteady flows behind a circular cylinder at conditions both near 
to and far from resonance. 
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3. Experimental setup 

The experiments were carried out in the multi-discipHnary low speed wind tunnel of the 
University of Liege. The Aeronautical test section of the wind tunnel was used, which 
measures 2m x 1.5m x 5m (width x height x length) and has a maximum airspeed of 60m/s. 




Fig. 1. Laser sheet illuminating 2D section of flow around a cylinder 

A circular cylinder of 36mm diameter and 1.32m span was placed in the wind tunnel, 
supported by its mid-span point near the middle of the test section. The cylinder was made of 
aluminium tube and painted matt black. The cylinder is rigidly supported but is flexible itself. 
Its first symmetric bending mode has a frequency of around 70Hz. Therefore, it is expected 
that when the frequency of the Von Karman vortex street behind the cylinder matches the 
first bending frequency, the free ends of the cylinder will oscillate quite visibly. Away from 
resonance, the cylinder will be static. This setup is ideal for the purposes of the present 
investigation, as it allows the examination of the unsteady flow behind both a static and a 
vibrating object. 

4. PIV system setup 

The PIV system used for these experiments consists of the following components: 

1. A Litron LDY301-PIV Q-switched laser system. It is a dual power, dual cavity laser with 
a wavelength of 527nm, switching at lOOOHz. The two laser beams contain 2 x lOmJ of 
energy. 

2. Optical modules for producing a laser sheet. 

3. A Phantom V9.1 camera with a maximum resolution of 1600 x 1200 pixels at a frequency 
of IKHz and 6GB of internal memory buffer. 
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4. A timer box for synchronizing the laser with the camera. 

5. A seeding generator with 3 bar back pressure suitable for PIV particle generation. 

6. Dantec Dynamics Studio PIV data acquisition and analysis software. 

7. A high specification Personal Computer. 

In practice, the system is capable of capturing PIV data in a window of around 13 x 11cm at a 

frequency of lOOOHz. 

The laser sheet was placed on the side of the cylinder nearest to the working section's 

observation window and aligned with the airflow, so as to illuminate a 2D section of flow 

around the cylinder. The laser sheet position can be seen in figure 1. Notice that the centre of 

the laser sheet lies aft of the cylinder. The very sharp shadow under the cylinder is also worth 

noting. 

4.1 PIV results 

Figure 2 shows a snapshot of the illuminated particles around the cylinder. It can be seen 
that the laser illuminates the aft upper section of the cylinder itself (white arc) as well as 
seeded particles on the upper surface of the cylinder and in the wake. Evidently, there are no 
illuminated particles in the shaded area under the cylinder. 




Fig. 2. Illuminated particles around the cylinder 

The PIV system takes two such snapshots at a very short time interval, typically 1 — lOOOf/s. 
A region of interest (ROI) is defined in the snapshots. This ROI is further divided into 
subregions (e.g. 64 x 64 pixels). The motion of the particles inside each subregion of the first 
snapshot is correlated to the second snapshot. The aim of the analysis is to decide where each 
subregion of photo 1 has moved to on photo 2. Thus, a velocity vector is placed in the centre 
of each subregion, representing the global motion of the particles inside the latter. The entire 
process is carried out by means of the Dynamics Studio PIV software and requires calibration 
data consisting of the free stream airspeed or a characteristic length, in this case the cylinder 
diameter. 

The end result of the PIV data reduction process is a velocity vector field calculated at each 
instance in time for which the visualization took place, as shown in figure 3. The image 
correlation process can sometimes lead to the calculation of bad vectors (outliers); these are 
detected and replaced by averages of all the neighboring vectors. 

PIV visualizations for the circular cylinder were carried out at airspeeds from 10 to 20m /s, at 
a sampling frequency of lOOOHz and sampling times from 0.1s to 4s. The recovered unsteady 
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Fig. 3. Computed velocity vectors for a sample snapshot 

vector fields had a resolution of 24 gridpoints in the y-direction and 41 in the x-direction. 
Therefore, the sizes of the u and v matrices ranged from 24 x 41 x 100 to 24 x 41 x 4000. 

4.1.1 Flow frequency variation with airspeed 

In order to verify that the PIV system and POD decomposition analysis are performing 
correctly, a large number of PIV measurements were carried out at airspeeds between 4m /s 
and 26m /s. The flow frequencies recovered by the POD method for all these measurements 
were then compared to the theoretical frequencies, assuming that the cylinder has a Strouhal 
number of 0.2. This comparison can be seen in figure 4, where the experimentally estimated 
frequencies are plotted as stars with error bars and the theoretical frequency is plotted as a 
dashed sloped line. The error bars represent the frequency increment, equal to the sampling 
frequency divided by the number of time measurements M and is equal to 2.02Hz. 
Figure 4 shows that the frequencies estimated from the decomposed PIV measurements 
are in good agreement with the theoretical predictions. It can be concluded that both the 
instrumentation and the POD analysis were correctly operated. The two dotted lines in the 
figure represent the airspeed limits at which significant cylinder vibration amplitudes were 
observed. Indeed, resonance phenomena were observed at airspeeds between 13.5m/s and 
18m/s, corresponding to vortex shedding frequencies of 70-105Hz. The lock-in phenomenon, 
whereby the flow frequency adapts itself to the structure's natural frequency throughout the 
resonance airspeed range, is not evident in this data. The reason for this absence of lock-in 
is that the measurements used for constructing figure 4 were taken close to the cylinder's 
midpoint, as seen in figure 1. At this location the amplitudes of the vertical vibration are 
small and have no impact on the shedding process. The frequency of the latter follows the 
linear Strouhal relation. 



4.2 POD analysis 

At the end of the PIV data treatment, a set of u (x, y, t) and v{x, y, t) matrices were obtained for 
each tested airspeed. These matrices where then analyzed using the POD procedure described 
in section 2. Sample results from three airspeeds are presented and discussed in this section. 
These are labeled as: 

• Test 1: Free stream airspeed of 18.8m/s, sampling frequency of lOOOHz, sampling time of 
0.298s, PIV laser sheet at 0.2m from the cylinder's midpoint. 
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Fig. 4. Comparison of estimated frequencies to theoretical frequencies 
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Fig. 5. Mean flow vectors for test 1 
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Fig. 6. Eigenvalue ratios for the first 20 eigenvalues for test 1 

• Test 2: Free stream airspeed of 13m/s, sampling frequency of lOOOHz, sampling time of 
0.099s, PIV laser sheet at 0.2m from the cylinder's midpoint. 

• Test 3: Free stream airspeed of 13.9m/s, sampling frequency of lOOOHz, sampling time of 
0.099s, PIV laser sheet at 0.2m from the cylinder's midpoint. 

• Test 4: Free stream airspeed of 14.8m/s, sampling frequency of lOOOHz, sampling time of 
0.099s, PIV laser sheet at 0.4m from the cylinder's midpoint 

The Reynolds numbers for the three tests range from 40,000 to 55,000, which puts the flowfield 
in the 'transitional in the shear layer' category. In other words, the boundary layer on the 
cylinder's surface is expected to be laminar. After separation, there is periodic ejection of 
vortices in the wake, as in the case of laminar vortex shedding, but the shear layer causing 
this ejection is transitional, giving rise to small turbulent eddies. 

Test 1 is used as the reference test, as it lies very far from aero-structural resonance and, 
therefore, there is negligible cylinder movement. Test 4 lies right on resonance and the 
cylinder vibrates significantly at the PIV measurement position. For all the tests, the first 
step in the POD procedure was to define the region of interest, so as not to include in the POD 
calculations the velocity vectors under the cylinder, which are not observable. The region of 
interest was therefore limited to a point just downstream of the cylinder. 

4.2.1 Test 1 

The next step in the POD procedure for Test 1 was to calculated the mean flow. This calculation 
involved the time averaging of the u{x,y, t) and v{x,y, t) matrices, leading to the mean flow 
shown in figure 5. It can be seen that the mean flow is essentially an area of slow, recirculating 
flow, located just behind the cylinder. In other words, it can be seen as the steady wake, to 
which an unsteady wake is superimposed. 

Once the mean flow was subtracted from u{x,y,t) and v{x,y,t) the unsteady vector fields 
u'{x,y,t) and v' {x,y,t) were obtained and Proper Orthogonal Decomposition was applied. 
The eigenvalue ratios A^ / T^^i K obtained for the first 20 modes are shown in figure 6. It can 
be seen that only the first two eigenvalues have significant contributions to the total energy in 
the unsteady flow, of 35% and 25% respectively. All higher eigenvalues have contributions of 
less than 5% and can therefore be neglected. 

The mode shapes associated to the two retained eigenvalues are shown in figure 7, plotted as 
filled contour plots. Subfigure 7(a) plots the values of (pu{x,y) (left) and (py{x,y) (right) for 
mode 1 as a filled contour plot, black signifying a low value and white signifying a high 
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Fig. 7. The first two mode shapes for test 1 

value. Subfigure 7(b) depicts the same information for mode 2. The horizontal distance 
between a maximum and a minimum grows from 30mm to 40mm with distance downstream. 
The vertical distance increases with downstream distance from 30mm to 46mm. If the mode 
shapes are assumed to be periodic, i.e. repeatable further downstream, then the two modes 
are separated by 1/4 of a period. It should be noted that the distances between maxima and 
minima are not constant because the measurement is just behind the cylinder, where there are 
large variations in the mean flow, as seen in figure 5. It is likely that the distances between 
maxima and minima are stabilized further downstream. 

Figure 8 shows the variation in time of the two retained generalized coordinates. It 
can be seen that the two generalized coordinates have the same fundamental frequency 
of 107.5Hz. Furthermore, they are both subjected to a beating phenomenon, with the 
response amplitudes dropping momentarily at around 0.12s and again at around 0.3s. This 
beating demonstrates that the flow is quasi-periodic, with significant variations in amplitude 
occurring momentarily. This quasi-periodic nature is justified by the fact that the flow is 
transitional in the shear layer and there are small turbulent eddies absorbing some of the 
flow energy. 
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Fig. 8. Variation in time of the two retained generalized coordinates. Test 1 

In fact, the two mode shapes of figure 7 are the dominant modes and can be viewed as 
Taminar' modes. The modes that have been neglected can be viewed as 'turbulent' modes, 
which contain little energy over the complete time history but can momentarily absorb energy 
from the dominant modes. This is exactly what happens in the case of this test. At time indices 
of 0.12s and 0.3s, the response amplitude of modes 1 and 2 drops significantly; simultaneously 
the response amplitude of mode 3 increases visibly, as shown in figure 9(a), which plots the 
variation of generalized coordinate cj^{t) with time. It can be clearly seen that the maxima of 
this mode occur at 0.12s and 0.3s. 

The mode shape for mode 3 can be seen in figure 9(b). It is clear that this mode shape is 
significantly noisier than the mode shapes of the first two modes, shown in figure 7. This 
noisiness is consistent with the hypothesis that mode 3 represents turbulent flow energy. It 
should be mentioned that higher modes do not demonstrate any clear increases in amplitude 
at 0.12s and 0.3s. Therefore, a complete model of the flow of Test 1 will contain: 

• two modes if only the laminar flow is of interest 

• three modes if it is desired to account for some of the energy momentarily lost by the 
laminar modes 



4.2.2 Test 2 

Test 2 is similar to Test 1 in the sense that resonance is not occurring yet; however, the 
condition is much closer to resonance than Test 1 and the amplitude of vibration of the cylinder 
is small but noticeable. Under these circumstances, the effect of evaluating a mean flow and 
subtracting it from the total vector field must be revised. In effect, as the mean flow is the 
steady wake behind the cylinder, if the cylinder oscillates by a significant amount, then the 
concept of a steady wake is no longer valid and the mean flow should not be evaluated. 
In other words, while it is still possible to calculate values for u and v over the entire time 
history, these values will not be the same over a different time history. In essence, the 
movement of the cylinder is an excitation force that is applied to the fluid; it will respond 
at the excitation frequency and at higher harmonics but there will be no response component 
with zero frequency. 

Nevertheless, if the amplitude of oscillation is very small, then evaluating and subtracting the 
mean flow will not cause big errors in the POD procedure. For Test 2, the POD method was 
applied twice, the first time after subtracting the mean flow and the second after subtracting 
only the wind tunnel free stream, U. In other words, in the first application the POD analysis 
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Fig. 9. The third mode for test 1 



was carried out on u^ and v^ using equations 3 to 7 while in the second it was carried out on 
u — U and v, such that 



^\hrh) = J^ J J {{u{x,y,ti)-U){u{x,y,t2) -U) + v{x,y,ti)v{x,y,t2))dxdy (8) 



with 



M 



PuA^^y) 



V^i^ m=l 



Y,{u{x,y,tni) -U)a^ 



PvA^ry) 



M 



/A^™^i 



^ v{x,y,tm)aM,i 



and 



'?<(0 = j J {{u(x,y,t) - U)(pu,i(x,y) + v(x,y,t)(p^^i{x,y)) dxdy 



(9) 



(10) 



For the first appHcation, the number of retained modes was 2, i.e. only the laminar unsteady 
flow modes were found to be significant; the resulting mode shapes were qualitatively 
similar to those obtained for Test 1 (figure 7). For the second application, the number of 
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Fig. 10. The first mode shape for test 2 (top) compared to the mean flow (bottom) 

retained modes was 3. The first mode represented the mean flow while the other two 
modes represented the laminar unsteady flow and were very similar to the modes of the 
first application and, consequently, of figure 7. It is interesting to compare the mean flow 
subtracted from the data in the first application to the first mode obtained from the second 
application. Figure 10 shows contour plots of the first mode (p^^i and (pu^2 evaluated from the 
application of POD to u — U and v (top two plots) and of the mean flow components U and 
V (bottom plots). It can be seen that the two sets of contour plots are very similar. Therefore, 
the POD procedure described by equations 8 to 10 will calculate the mean flow as the most 
energetic mode. 

This is quite an interesting result because it suggests that there is no need to subtract the mean 
flow. If there is a steady component to the flow, then it will be identified automatically as the 
first mode. If there is no steady component, then equations 8 to 10 should be used anyway. 
Figure 11 shows the time response of the three resulting generalized coordinates. It can be seen 
that the generalized coordinate of the first mode is nearly constant with a value of around 600, 
while the coordinates of the other two modes oscillate around zero. It is clear that the POD 
procedure can differentiate between steady and unsteady responses. 
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Fig. 11. Variation in time of the three retained generalized coordinates. Test 2 

4.2.3 Test 3 

For Test 3 the POD procedure was carried out only using equations 8 to 10. As mentioned 
before. Test 3 gave rise to significant amplitudes of cylinder vibration as the flow frequency 
was very close to a natural frequency of the cylinder. Again, three modes were retained, one 
representing a mean flow and two representing the laminar oscillating flow. 
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Fig. 12. Variation in time of the three retained generalized coordinates. Test 3 

The response amplitudes of the three generalized coordinates were significantly higher than 
in the case of Test 2, as seen in figure 12. Even the generalized coordinate of mode 1 is far more 
unsteady, although its variation in time is not periodic. This aperiodic variation suggests that 
there is a component of the cylinder vibration in the flow data but it is of the same order 
as the turbulent disturbances and /or experimental error. This was indeed the case, as the 
measurement point was close to the cylinder's midpoint (see figure 1), therefore the local 
vibration amplitude was small, of the order of less than 1mm. 



4.2.4 Test 4 

The final test was carried out at a slightly higher airspeed but, crucially, the PIV laser sheet 
was positioned at a span- wise point further from the cylinder's midpoint than for the other 
three tests. At this particular span- wise position, the cylinder was vibrating significantly, with 
an amplitude of nearly 2mm and a frequency of 70.5Hz. 

The POD procedure was again applied using equations 8 to 10, i.e. without subtracting the 
mean flow. The resulting mode shapes were similar to those obtained during Test 3. The 
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Fig. 13. Variation in time of the three retained generalized coordinates. Test 4 

resulting generalized coordinates are plotted on figure 13. It can be clearly seen that the 
response of the first generalized coordinate is now much more oscillatory than in the case 
of figure 12. Furthermore, despite the randomness of this response, there is a clear periodic 
component at a frequency close to that of the oscillating modes. 
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Fig. 14. Power Spectral Density plots of the generalized coordinates and the cylinder 
displacement. 

The observations of Test 4 suggest that POD can decompose flowfields that feature 
unsteadiness due to the presence of both separated flow and structural motion of the wind 
tunnel model. However, for this decomposition to be successful, the response amplitude 
of the structure must be significantly higher than the size of the turbulent eddies. In such 
cases, straightforward POD will result in generalized coordinates that contain the response 
frequencies of both the separated flow and the structural motion. Figure 14 shows Power 
Spectral Densities (PSD) of the qi{t) and z{t) signals, z{t) being the vertical displacement 
time history of the cylinder at the PIV measurement position. The PSDs were calculated 
using the Welch method, with a Hamming window 512 samples long and 50% overlap. The 
cylinder's displacement response clearly contains only one frequency component at 70.5Hz. 
The generalized coordinates feature two frequency components, one at 70.5Hz and a stronger 
one at 85.9Hz. It can be inferred that 70.5Hz is the structural response frequency while 85.9Hz 
is the vortex shedding frequency. The Strouhal frequency at the Test 4 airspeed is 82.2Hz if a 
Strouhal number of 0.2 is assumed (see figure 4), i.e. quite close to 85.9Hz. 
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The generalized coordinates can be seen as the responses in time of the fluid due to both 
flow unsteadiness and cylinder motion. Furthermore, the cylinder motion can be seen as 
an external excitation acting on the fluid. Therefore, it possible to set up an input-output 
POD model, whereby the input is the cylinder motion and the outputs are the generalized 
coordinates. Frequency Response Functions (FRF) can then be created of the form 



Hi(co) 



Qijco) 
Z{co) 



(11) 



where Hj is the zth FRF, Qi is the fth generalized coordinate in the frequency domain, Z 
is the cylinder displacement in the frequency domain and co is the radial frequency. Such 
FRFs can be estimated using a Welch-type windowed approach and involving cross and 
auto-correlations of the outputs and input. 
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Fig. 15. FRFs between the generalized coordinates and the cylinder displacement. 

Figure 15 shows the FRFs estimated for the first three modes of Test 4. It can be seen that 
the main frequency component of all FRFs is the vortex shedding frequency at 85.9Hz, in the 
frequency range between 50 and lOOHz. The FRF of mode 1, Hi, features several other peaks 
but these are most likely noise; for this mode the main frequency component lies at OHz, since 
the mode mainly reflects the mean flow. 



4.2.5 Mode shapes variation with airspeed 

Here, the variation of the modes shapes with airspeed will be discussed. Only the free stream 
is subtracted from the local velocity field for the results in this subsection, not the mean flow. 
It was observed that the first three modes, i.e. the steady flow and the two laminar modes, 
change very little in shape with airspeed. 
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(a) Mode 1 



(b) Mode 2 




(c) Mode 3 



Fig. 16. Vector plots of the first three mode shapes 

Figure 16 shows plots of the first three modes in vector form, whereby the vectors' horizontal 
component is taken as (puj{x,y) and the vertical one as (pyj{x,y). Mode one (figure 16(a)) 
consists of two areas of recirculation behind the top and the bottom of the cylinder. Between 
them, they cause a significant area of flow towards the cylinder with a height approximately 
equal to the cylinder's diameter. Mode two (figure 16(b)) consists of a large vortex, positioned 
approximately one diameter behind the cylinder and centered on the cylinder's centerline. 
It is accompanied by two much smaller areas of recirculation lying just behind the top and 
bottom of the cylinder. Mode 3 (figure 16(c)) consists of two counter-rotating vortices at 
approximately half a diameter (the strongest) and one diameter (the weakest) behind the 
cylinder. They are both lying on the cylinder's centerline. As mentioned above, these three 
modes remain largely unaffected by the airspeed. 

Figure 17 shows vector plots of mode 4 obtained for the four different tests, i.e. at four different 
airspeeds. It can be seen that all mode shapes are different. The mode shape for Test Iconsists 
of four vortices; for Test 2 of three vortices, one of which is very weak; for Test 3 the mode 
shape consists of three vortices again but in a different arrangement; finally, for Test 4 the 
number of vortices is difficult to determine because they are quite weak and small. 
The reason for the change in mode shape 4 with test airspeed may be numerical rather than 
physical. In fact, mode shapes very similar to that of figure 17(a) appeared as mode 6 in 
Test 3 and mode 5 in Test 4. Therefore, the same mode shape can be preserved but its 
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Fig. 17. Vector plots of the fourth mode shape for the four different tests 

relative importance changes between tests. This phenomenon may provide justification for the 
assertion presented here concerning the higher modes, i.e. that modes higher than 3 represent 
transitional effects. 



5. Conclusions 

The feasibility of applying Proper Orthogonal Decomposition to experimentally measured 
flows around vibrating structures has been demonstrated. It has been shown that this type 
of decomposition analysis can provide some very interesting data about the observed flows, 
such as the dominant mode shapes and frequencies. Furthermore, it was shown that structural 
vibrations can be detected by the POD procedure applied on PIV flow visualization data using 
an output-only approach. 

By considering the cylinder structural response as a forcing function, it is possible to 
create input-output POD models, whereby the generalized coordinates can be obtained from 
Frequency Response Functions relating the cylinder displacement response to the generalized 
coordinates themselves. It is shown that such FRFs feature two main frequency components, 
the mean flow frequency (i.e. OHz) and the vortex shedding frequency. Therefore, they are 
independent of the structural response frequency. 
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1. Introduction 

Control surfaces, such as trailing-edge flaps, provide a means to dynamically alter the 
aerodynamic characteristics of aircraft for primary flight control, secondary vibration 
control, and even higher frequency noise control. While the development of several novel 
technologies has been explored, many practical implementation barriers still exist for a 
single actuation system to serve all three of these objectives. This is particularly true for 
rotor craft, where the demands of the harsh rotary and vibratory environment are severe in 
terms of actuator force and displacement, bandwidth limitations, life cycle concerns, and 
physical volume available. Accordingly, it has been assumed that on-blade active controls of 
a rotorcraft would be subject to the most stringent requirements in the subsonic flight 
regime, and if a control surface actuation technology could survive here, it could be 
reasonably applied to a fixed wing aircraft. A brief account of the current state-of-the-art for 
rotorcraft blade controls follows. 

Helicopter rotors typically operate in a highly unsteady aerodynamic environment. In forward 
flight, the rotor blade sections experience large variations in angle-of-attack over one 
revolution. This is the primary source of performance degradation, such as high vibration and 
retreating blade stall. Actively changing the angle-of-attack of the blade sections as a function 
of blade azimuth has been shown to significantly alleviate vibration levels, as well as improve 
aerodynamic performance of the rotor (Straub et al., 2000). The change in angle-of-attack can 
be accomplished in a variety of ways. Implementation of high bandwidth hydraulic actuators 
in the rotating frame has demonstrated the ability to actively change the root pitch of the rotor 
blades and has since been demonstrated in both scale models (Lorber et al., 2001) and full scale 
tests (Arnold & Strecker, 2002). Another approach is to vary the aerodynamic forces on the 
blades by dynamically changing the geometry of the airfoil sections. This can be accomplished 
through actively controlling blade twist, airfoil camber, or through the use of trailing-edge 
flaps (Hall & Wereley, 1993). Recent advances in adaptive materials have led to a variety of 
schemes for on-blade actuation in these areas (Chopra, 2000). Some of these include 
piezoelectric innovations such as adaptive twist of the rotor blade (Chen & Chopra, 1997; Chen 
et al, 2001; Shin et al, 2005), trailing-edge flaps (Straub et al, 2001; Fulton, 2000; Fulton, 2005), 
and active camber control (Konstanzer et al., 2001; Nissly et al., 2005). 

In comparing these approaches to active rotor systems, there are potential drawbacks, 
however. For instance, implementation of hydraulic systems in the rotating frame of 
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production helicopters is a challenging task due to the complexity of the system, the 
increase in maintenance associated with the large number of moving parts, as well as the 
hydraulic slip ring, and the added on-blade mass associated with the weight of the 
hydraulic fluid and piping. Considering the active blade twist concept, there is also a large 
weight penalty due to the distributed nature of the actuators. There are also several 
unanswered questions as to whether active material solutions like piezoelectrics can survive 
in the operational environment, whether they have sufficient fatigue life for practical 
consideration, or whether they are properly scaled for operating over broad deflection and 
frequency ranges in full size rotors. Despite this, there have been numerous developments 
aimed at demonstrating the potential that active materials have for full scale rotors. Most 
notable is the piezoelectric actuator work conducted independently by Boeing and 
Eurocopter. Both using piezoelectric actuators to drive trailing-edge flaps on rotor blades, 
the Boeing development led to full scale whirl testing for vibration control (Straub et al., 
2004) and a full scale wind tunnel test for noise reduction (Straub et al., 2009), and the 
Eurocopter development led to full scale flight testing for vibration control (Roth et al., 2006; 
Dieterichetal.,2006). 

Trailing-edge flaps provide localized actuation and can generate significant control 
authority, but these discrete control surfaces do increase drag from the discontinuities and 
abrupt changes in airfoil contour. The active camber control concept alleviates this issue by 
varying the camber of the airfoil to produce a conformal shape change. There are several 
technical barriers that exist in actual implementation, such as the development of a flexible 
skin and supporting core that can withstand the harsh rotating environment of a helicopter. 
These topics are beginning to be addressed in the fixed frame, such as in unmanned air 
vehicles (Flanagan et al, 2007; Bubert et al, 2010; Olympio & Gandhi, 2010), but the 
technology has not reached the maturity level required for rotorcraft. Therefore, it appears 
that the trailing-edge flap may be the leading candidate control scheme for active rotors 
given the current state-of-the-art in practical actuation strategies. 

Using these drawbacks as motivation to investigate alternative approaches to active 
aerodynamic control, a less conventional, yet properly scalable, approach to trailing-edge 
flap actuation has been developed and tested, and it employs Pneumatic Artificial Muscles, 
or PAMs, as actuators (Kothera et al., 2010). These actuators, originally developed for 
orthotic devices in the 1950s by J.L. McKibben (Gaylord, 1958; Schulte, 1961), typically 
contract in response to an increase in internal pressure, and have been used in robotic 
applications (Tondu et al, 1994; Medrano-Cerda et al, 1995; Daerden & Lefeber, 2002). 
Typically constructed of an elastomeric bladder surrounded by a braided sleeve, the stroke 
in these low cost and light weight actuators comes from re-orientation of the stiff braid 
fibers as the bladder expands radially. Previous research with these devices has 
experimentally confirmed their applicability to trailing-edge flaps in both the fixed frame 
(Woods et al, 2007; Kothera et al, 2008; Woods et al, 2010a) and the rotating frame (Bubert 
et al., 2007; Woods et al., 2010b). The development of and results from these fixed-frame 
tests will be presented here from two different mechanical configurations and loading 
conditions, as a step toward demonstrating the feasibility of using PAM actuators for 
aerospace applications. The next section will discuss the general system design, which will 
be followed by a discussion of bench- top testing results. Then the wind tunnel test article 
development is shown, along with experimental test data. 
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2. System design 

Three different PAM actuation systems were designed, built, and tested, with experimental 
evaluations taking place first on the bench and then moving on to a wind tunnel. In each 
case, the aerodynamic hinge moment was predicted first and then the design of the 
actuation mechanism followed. 

2.1 Aerodynamic prediction 

The actuation system design used two-dimensional thin airfoil theory to predict the hinge 
moments that were used in sizing the actuators and bench-test loading springs to simulate 
the aerodynamic hinge moments. Environmental parameters used in the calculations are 
those of Sea Level Standard (SLS), which have a temperature of T = 288 K, pressure of P = 
101.3 kPa, and air density oi p = 1.225 kg/m^. The maximum speed of the Glenn L. Martin 
wind tunnel (GLMWT) at the University of Maryland is Mach 0.3, so this was also used in 
the design as a maximum free-stream condition, and a reasonable angle-of-attack of a = 6^ 
was also incorporated. It was assumed in this loading mechanism design that there was no 
flow separation over the airfoil and that the gap between the flap and the airfoil was sealed. 
The most important quantity for prediction in terms of sizing the actuation system is the 
hinge moment. This can be calculated according to 

H = -pV^SfCfCj^, 

where Ch is the hinge moment coefficient, Cha is the hinge coefficient due to angle-of-attack, Chf 
is the hinge coefficient due to flap deflection, H is the actual hinge moment, v is the wind 
speed, Sf is the flap area, and c/ is the chord of the flap. For the angle-of-attack effects, 
tabulated airfoil data can be used with the equation 



^hc 



= (J'i)o+2[(fli);^-(«i);]tan(0.5r-f/c), 



where (hi) is the rate of change of hinge moment coefficient with incidence, the (ai) terms are 
related to the lift-curve slope, x is the trailing-edge angle, and t/c is the thickness ratio (Etkin, 
1982). The effect on the hinge moment due to flap deflection is computed by 



dci, dc 



dcj dS 



^i 



where dcj/dci is the rate of change of the hinge moment coefficient with respect to the change 
in lift coefficient, dcj/d5 is the rate of change of the hinge moment coefficient with respect to 
the change in flap deflection, q is the lift coefficient, and 5 is the flap deflection (Abbott & 
von Doenhoff, 1959). Because this term depends on the lift coefficient, its computation was 
also required in the prediction. For the current case with a flap on the trailing-edge, added 
lifting effects from the flap on the lift coefficient can be approximated as 
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where k is the fraction for how the flap deflection angle changes the effective angle-of-attack 
of the airfoil (Eastman & Pinkerton, 1930). Here, the lift coefficient due to angle-of-attack is 
estimated with Mach scaling according to 

2;r 



where M is the Mach number. The denominator in this equation is commonly denoted as 
the parameter p. 

As was mentioned, a total of three PAM-flap systems were evaluated. Table 1 lists the key 
design parameters for each and defines each with a ''system number,'' which will be used 
throughout to more easily identify the results being shown. As can be seen, all three systems 
were sized for NACA 0012 airfoils with nominal angles-of-attack of 6 degrees. It should also 
be noted that each successive system was designed and tested for a higher Mach number. 
The first two systems were identical in terms of the overall model size and flap size, with the 
second containing stronger PAM actuators for the higher air loads. The third system was 
designed for an even higher wind speed and a smaller airfoil section to better illustrate the 
capability of the PAM actuators for scalability and high performance. 



System No. 


1 


2 


3 


Wind Speed (M) 


0.1 


0.3 


0.56 


Angle-of- Attack (deg) 


6 


6 


6 


Airfoil Type 


NACA 0012 


NACA 0012 


NACA 0012 


Airfoil Chord (in) 


21 


21 


10.5 


Flap Chord (in) 


3.15 


3.15 


1.61 


Flap Span (in) 


10 


10 


34 



Table 1. Details of specific trailing-edge flap systems 

Using the aerodynamic equations listed above with the airfoil specifics and wind speeds in 
Table 1 led to predictions of the hinge moments that would need to be generated by the 
PAM actuation systems to deflect the various flaps. Figure 1 displays the simulation results 
for the noted systems and conditions. In each set of results, the estimated, bi-directional 
hinge moment contributions are shown from the airfoil angle-of-attack (green) and the flap 
deflection (red), along with the summed total (blue). Figure 1(a) shows the hinge moment 
for system 1 at the design condition. Figure 1(b) shows the hinge moment for system 2 at the 
design condition, and Figure 1(c) shows the hinge moment for system 3 at the design 
condition. Figure 1(d) has also been included to estimate the wind tunnel test condition for 
system 3. While designed for Mach 0.56, the test facility has a maximum wind speed of only 
Mach 0.3, so the designed actuation system was bench tested to the full spring-simulated 
aerodynamic stiffness, but it was wind tunnel tested to a reduced load with input pressure 
scaled accordingly. In order to meet the goal of at least 10 degrees of flap deflection, the 
PAM actuators must be able to produce 2.2 in-lb, 23 in-lb, and 70 in-lb for the noted design 
conditions of systems 1-3, respectively. Being that system 3 would be designed for 
operating at Mach 0.56, it would have little trouble meeting the reduced hinge moment 
requirement of only 22 in-lb for the wind tunnel test. 
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Fig. 1. Hinge moment predictions for actuation system designs - (a) system 1; (b) system 2; 
(c) system 3 bench-top; (d) system 3 wind tunnel 

2.2 Actuator design 

For the most part, the PAM actuation systems were sized to fit entirely within the airfoil 
contour. A schematic diagram of how the PAM system functions is provided in Figure 2. In 
this figure, the upper PAM is inactive (open to atmosphere. Pi = 0) and the lower PAM is 
active (pressurized, P2 > 0). The basic operation of a PAM actuator is to produce a contractile 
stroke upon internal pressurization. This stroke is the result of radial expansion of the 
elastomeric bladder, which is surrounded by a double helical weave of stiff fibers (e.g., 
polyethylene terephthalate - PET). As the softer bladder expands, the stiff fibers maintain 
their length and reorient themselves to accommodate the radial expansion. Consequently, 
the length of the device decreases. Also as indicated in the figure, any time pressure in one 
of the PAMs exceeds that in the other antagonistic PAM, a moment is generated about the 
flap hinge. The case shown is Pi < P2, which causes the flap to deflect in the downward 
direction. 

There are two key equations for predicting the deflections of the antagonistic PAM actuation 
system. The first is the arc length formula 



AL L^ - L2 



-S, 
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where Li is the length of the inactive PAM, Lz is the length of the active PAM, and r is the 
hinge radius. These are all labeled in Figure 2 with the convention that downward flap 
deflection angles are positive and upward flap deflection angles are negative. 




Fig. 2. Diagram of bi-directional PAM-flap actuator 

The second equation considers the available actuation force by 



AF = R 



r 



where Fi is the inactive PAM force (e.g., at psi), F2 is the active PAM force, and H is the 
hinge moment generated on the flap. For bi-directional operation of the flap, an antagonistic 
actuator arrangement is typically preferred. 

Based on these actuator design equations and the aerodynamic hinge moment predictions 
from section 2.1, an in-house database of PAM actuator performance data was consulted to 
select the actuator characteristics best suited to each of the three systems. Table 2 displays 
the results of this actuator sizing analysis. For the first two systems, which are essentially the 
same except for the increased strength of the second, a chordwise PAM orientation was 
selected. This means that the length and stroke direction of the PAM actuators is along the 
airfoil chord, which was made possible by the large airfoil section (21-in chord). System 3 
could not use chordwise actuators because of the higher forces required and the smaller size 
of the airfoil section. In this case, a spanwise orientation was selected, whereby the spanwise 
pulling motion of the PAM actuators is transferred into chordwise motion to deflect the flap 
through the use of a kinematic mechanism. While adding complexity and increasing the 
part count, having a mechanism does allow for consideration of mechanical advantage 
trade-offs to better tune system performance. For instance, when comparing the system 2 
and 3 actuators, it can be seen that larger diameter PAMs were used for system 2 although 
the expected loads were smaller here than for system 3. Note that it is a general design rule 
for PAMs with the same braid angle that larger diameter actuators will produce more force 
for a given pressure. This use of smaller diameter PAMs of essentially the same braid angle 
was facilitated by the mechanism, where stroke could be traded for force. As indicated in 
the table, the PAMs for system 3 are over 9 inches long, whereas the PAMs for system 2 are 
just over 3 inches in active length. 

Force-contraction data was collected from experiments conducted on a servo-hydraulic load 
frame. For each case, the PAM actuator was installed in the machine at its resting length, as 
shown in Figure 3(a). Next, the noted pressure was applied and held constant throughout 
the test. Immediately after the pressure was applied, the actuator blocked force 
measurement was recorded. Then, the actuator was allowed to contract slowly, or quasi- 
statically, to the point where no force was measured. This point is known as the free 
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System No. 


1 


2 


3 


Orientation 


Chordwise 


Chordwise 


Spanwise 


No. PAM Pairs 


1 


2 


1 


Diameter (in) 


0.5 


1.0 


0.625 


Active Length (in) 


3.15 


3.15 


9.10 


Resting Braid Angle (deg) 


47 


60 


61 


Braid Material 


PET 


PET 


Kevlar 


Bladder Material 


Latex 


Latex 


Latex 


Mechanical Advantage 


1 


1 


1.15 



Table 2. Details of specific PAM actuators selected for system designs 






(a) (b) 

Fig. 3. Characterization testing of PAM actuator - (a) resting length; (b) free contraction 

contraction value and is shown in Figute 3(b). After recording this value, the PAM actuator 
was stretched back to its resting length and the test was stopped. Typical performance data 
from each of the three selected PAM configurations is provided in Figure 4 at various 
pressure settings. Note that the x-axis data is presented as the non-dimensional contraction 
percentage. For the range of PAM lengths considered in this work, it has been shown that 
contraction percentage is largely independent of PAM length for a given diameter and braid 
angle. As an example, a particular PAM at a given pressure may have a free contraction of 
25%. This actuator with a 4-in resting length would thereby have a free contraction of 1 inch, 
whereas a 12-in resting length actuator would have a free contraction of 3 inches. 
It can be seen in Figure 4 that both force and contraction increase with pressure. For the 
performance metric of blocked force (measured force with 0% contraction, i.e., resting length 
and maximum possible force), the increase in force is linear with pressure, but the same is 
not true for the performance metric of free contraction (measured contraction with lb force, 
i.e., maximum possible stroke). The free contraction increases with pressure, but tends to 
roll off as a maximum limit is approached. These trends are valid, of course, only over the 
pressure range which does not lead to yield in the braid fibres. It should be noted here that 
our tests have shown that the PET braid does not begin to yield until approximately 200 psi 
and the yield point for the Kevlar braid is over 1000 psi. Given that the three system designs 
here are all intended to operate below 100 psi, a more than adequate safety margin exists for 
burst failure. Another characteristic to mention in regard to the typical force-contraction 
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Fig. 4. Performance data from PAM actuators selected for system designs - (a) system 1; (b) 

system 2; (c) system 3 

behaviour of PAM actuators is that some measure of hysteresis does exist in the loop. 
Despite the noticeable hysteresis, average values and a linear force-contraction 
approximation between blocked force and free contraction may be sufficient for initial 
system design and component sizing exercises. To state the maximum values of blocked 
force and free contraction (at 90 psi) for the system 1-3 actuators, respectively, we have 24 
lb with 14%, 250 lb with 27%, and 200 lb with 29%. 



3. Bench-top testing & validation 

3.1 Experimental setup 

To evaluate the PAM trailing-edge flap actuation systems prior to entering a wind tunnel 
environment, each of the three actuation systems was tested on a laboratory bench. Figure 5 
shows the test setups that were designed, fabricated, and tested. First in Figure 5(a) is a 
system sized for low subsonic air loads with the single pair of antagonistic PAM actuators 
(system 1) oriented along the chord of the airfoil. Here, the compressed air enters and exits 
the PAM actuators from their end near the upper-right corner of the photograph. Two 
aluminum extensions are at their other end and are instrumented with resistive strain gages 
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for measurement purposes. Representative inertia is mounted on the axis of flap rotation 
and a bending spring is attached to simulate aerodynamic stiffness. Figure 5(b) shows a 
similar arrangement with chord-oriented PAM actuators, but this system has two 
antagonistic pairs of larger PAM actuators (system 2) to increase the quasi-static flap 
deflection output. It can be seen in the figure that the PAMs attach to the hinge from one 
side and a spring loading mechanism attaches to the hinge axis from the other side. Linear 
extension springs were used for this test. Figure 5(c) shows an antagonistic pair of PAM 
actuators intended for spanwise orientation in the airfoil section, which uses a mechanism 
(not pictured) to turn the actuator work into chordwise motion that rotates the flap about its 
hinge axis (system 3). An inertial element was attached to the top of the axis of rotation here 
and linear extension springs were connected at a different radius than the PAMs to account 
for the mechanism that would be installed in the airfoil section. Also note that the difference 
in physical appearance of the PAMs from system 1 and 2 (black) to those in system 3 
(yellow) is due to the use of a different braid material. 









(a) (b) (c) 

Fig. 5. Experimental setups for bench-top evaluations - (a) system 1; (b) system 2; (c) system 3 

For each of the three test cases, the hinge rotation angle was measured with a Hall-effect angle 
sensor (Midori, series QPC). Two of the PAM actuators were instrumented with load cells 
(Honeywell) and the actuation pressure was measured with a pressure transducer (Omega, 
series PX209). Solenoid valves were used to direct the air flow into and out of the PAMs, as 
directed by a square wave voltage input for these open-loop experiments. For the two 
chordwise PAM configurations, the control valves used were SMC model VZ solenoids and 
for the spanwise PAM configuration, it was a Festo model MPYE valve with a much higher 
flow capacity. A National Instruments data acquisition system and laptop computer were used 
to run the experiments and to collect data. The test procedure included driving the control 
valves at various input frequencies and under different spring loads to measure the dynamic 
response of the PAM actuation systems in terms of flap deflection angle bandwidth. 



3.2 Bench-top test results 

Figure 6 shows the bench-top test results from the three cases as half peak-to-peak deflection 
values that were averaged over several actuation cycles. Here, Figure 6(a) shows data from 
system 1 for a range of pressures. It should be noted that the spring load indicated in the 
caption simulates air loads existing at a free-stream velocity of Mach 0.23, while the design 
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condition was only for Mach 0.1. The difference here was due to spring availability at the 
time of the test, so a more conservative loading condition (approximately four times the 
expected wind tunnel load) was tested on the bench- top. Given that the PAM actuation 
system was able to produce nearly ±5 degrees of flap deflection over the tested frequency 
range of 30 Hz provided sufficient assurance that the system would meet, and most likely 
surpass, the ±10 degree target in the wind tunnel. As expected from basic PAM operation, 
this figure also shows that increasing the pressure led to increased deflection capability. 
Figure 6(b) shows data from the Mach 0.3 spring load of the chordwise, double PAM pair 
actuators (system 2) at various pressures, showing a different characteristic than was shown 
for system 1. As can be seen, the quasi-static deflection output increases fairly linearly with 
actuation pressure, but there is a sharp roll-off in deflection performance at frequencies 
above 1 Hz that reduces the achievable deflection level at higher frequencies. Additionally, 
the measured deflection at high frequency appears to be independent of pressure. These 
response features indicate that air flow restrictions were present in the pneumatic supply 
system. This is a reasonable outcome when considering that the same pneumatic valves 
were used for systems 1 and 2, while the air volume required for operating system 2 was 
much larger than that for system 1. This was later improved upon in a revised bench- 
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Fig. 6. Bench-top test results - (a) system 1 with Mach 0.23 spring load; (b) system 2 with 
Mach 0.3 spring load; (c) system 3 with Mach 0.56 spring load 
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top test, but the improved pneumatics for system 2 was not tested in a wind tunnel, so those 
results will not be presented here (Kothera et al., 2008). In continued scaling of the PAM 
actuation system, system 3 did incorporate improved pneumatic elements. Figure 6(c) 
shows experimental flap deflection measurements from a Mach 0.56 spring load at various 
pressures for system 3, where it can be seen that the deflection bandwidth was significantly 
increased over that seen from system 2, with ±10 degrees now being produced up to nearly 
10 Hz at 90 psi. Based on these results, it was concluded that systems 1 and 3 would meet 
the goal of at least ±10 degrees of flap deflection in the dynamic response for their respective 
wind tunnel tests, but system 2 would likely fall short of the goal due to the flow rate 
restrictions in the pneumatic circuit. These tests also verified proper basic functionality of 
the actuation systems. 

4. Wind tunnel testing 

After successful bench-top evaluations, the trailing-edge flap actuation systems were placed in 
various wind tunnel test articles to assess their performance under actual aerodynamic loading 
conditions. This section will discuss the test article fabrication and experimental test results. 

4.1 Wind tunnel test article development 

The three PAM actuation systems were integrated into different test articles based on their 
originally specified performance conditions and measured performance from the bench-top 
evaluations. Systems 1 and 2 had the same geometry, as stated in Table 1, but different 
construction was used to better match the test conditions. That is, system 2 was designed with 
a stronger frame structure and skin than system 1. The frame for system 1 was constructed 
with 0.25-in aluminum and consisted of a spar running the span of the airfoil section (24-in), 
an airfoil-shaped end plate at both ends, and two intermediate support ribs that extended from 
the spar to the trailing-edge. The flap hinge was fabricated from a precision ground 0.25-in 
stainless steel rod, and the flap was able to rotate about the hinge with a bearing at each end. 
The skin was made with a wet lay-up process using two plies of pre-impregnated fibreglass, 
and formed the shape of the airfoil section with a Styrofoam mold. Additional reinforcement 
at the leading-edge was provided with a thin sandwich structure of fibreglass and foam, which 
formed a shell-like structure around the spar so that the solenoid valves and air tubing could 
be installed inside. The system 2 wind tunnel model had the same basic layout as the system 1 
model, but it was made to be stronger since it was to be tested at Mach 0.3. That is, it had two 
airfoil-shaped end plates and two mid-span support ribs, all constructed of 0.5-in aluminum, 
and a steel spar. The entire skin was also fabricated as sections of sandwich structure 
composites with foam core inside 2-ply fibreglass. There was a separate leading-edge, or D- 
spar section, and aft panels that were fastened to the frame structure. As with system 1, system 
2 also placed the solenoid control valves inside the D-spar cavity surrounded by the shell-like 
skin. Figure 7(a) shows the system 1 wind tunnel model installed in the test fixture between 
two acrylic vertical walls and positioned in front of the free-jet wind tunnel. Figure 7(b) shows 
the wind tunnel model for system 2, which was fixed at the floor of the Glenn L. Martin wind 
tunnel at the University of Maryland and cantilevered upward. There was also an elliptical top 
plate bolted to the top of the cantilevered model to promote more desirable air flow 
characteristics over the test section. The primary difference in appearance of these two models 
is that the system 1 model was left bare and the system 2 model was painted black. Note also 
that strips of aluminum tape were placed over the recessed fasteners. 
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The wind tunnel model for system 3 was fabricated with an entirely different approach than 
those for systems 1 and 2. In this case, instead of building the complete airfoil section from 
raw materials, the outboard section of a scrapped full-scale Bell 407 helicopter rotor blade 
was employed as the basis for the wind tunnel model. Since this rotor blade does not 
currently feature a trailing-edge flap, the blade was physically modified to accept a flap and 
the PAM actuation system. The spanwise PAM actuators and mechanism were mounted 
inside the D-spar of the blade and the aft section of the blade was removed. This section was 
first used to create a mold for the new trailing-edge flap so that the geometry of the blade 
could be matched exactly. Note that the production Bell 407 blade contains taper and twist, 
which would have otherwise been difficult to match in constructing the flap. Aluminum 
strips were bonded inside the remaining aft portion of the blade skin where the flap section 
had been removed in order to provide bearing surfaces and structural support to fasten the 
trailing-edge flap in place. The flap was supported at the inboard and outboard ends and at 
two intermediate locations. A hole was drilled in the aft wall of the spar and a slot was 
placed in the skin to attach a control rod from the mechanism inside the spar to a control 
horn on the trailing-edge flap. These parts are on the opposite side of blade section in Figure 
7(c), so they are not visible in the photograph; the control rod was located at the inboard 
edge of the flap, or near the floor in the figure. Recall that system 3 was the only system for 
which a component of the PAM actuation system violated the baseline airfoil profile. This 
was due largely to the smaller size of the airfoil section and the higher wind loading 
condition. Also note that the tip of the cantilevered blade section was laterally restrained 
with a cable to limit side-to-side motion of the blade during flap actuation. 




Fig. 7. Wind tunnel test articles - 
3, closed-circuit 



(a) system 1, free-jet; (b) system 2, closed-circuit; (c) system 
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4.2 Wind tunnel test results 

Test results in the form of averaged, half peak-to-peak deflection are displayed in Figure 8, 
following the same layout as the previous figures. Shown in each case are the measured 
results from the previously specified wind tunnel test conditions, with the effect of varying 
actuation pressure also displayed for systems 1 and 3. Reaching a maximum speed of Mach 
0.1, Figure 8(a) shows the open-loop dynamic response of system 1 at various actuation 
frequencies and pressures. As was inferred from the bench- top test, this PAM actuation 
system was able to far exceed the original goal of ±10 degrees at high frequency. At only 30 
psi operating pressure in the PAMs, this system was able to produce ±10 degrees beyond 20 
Hz, and operating the PAM actuators with 90 psi led to ±20 degrees of flap deflection being 
produced up to nearly 25 Hz. There is also a resonance phenomenon apparent in this data 
set, which can be seen to increase in frequency with pressure. This changing resonance 
frequency is attributed to the changing stiffness of the PAM actuators as their operational 
pressure changes. Figure 8(b) shows the experimental flap deflections from system 2 that were 
measured at 90 psi in the PAM actuators and at two different angles-of -attack, though there is 
little noticeable difference in the dynamic response of the system at the two angles-of-attack. 



■ 90 psi (6°) 
-90 psi (12° 




Fig. 8. Wind tunnel test results - (a) single PAM pair, chordwise at Mach 0.1; (b) double 
PAM pair, chordwise at Mach 0.3; (c) single PAM pair, spanwise at Mach 0.3 
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As was expected from the bench-top test results, this system again illustrates a rapid drop 
off in achievable flap deflection as the actuation frequency increases. Recall that this was 
due to flow limitations in the pneumatic components. The ability to produce almost ±40 
degrees quasi-statically at Mach 0.3, however, is a promising result for the technology, 
especially when the dynamic response shown can be viewed as potentially a worst case 
situation achieving ±4 degrees of flap deflection up to 40 Hz. 

Figure 8(c) shows the wind tunnel results for system 3 at Mach 0.3. Recall that this is 
reduced from the bench-top test condition (Mach 0.56), but was the maximum possible 
speed of the wind tunnel used for testing. There are also two lines for each of the noted 
actuation pressure levels. The solid line represents the flap deflection measured at the 
inboard edge of the flap and the dotted line represents the deflection at the outboard edge of 
the flap. Since there is some difference between the two ends of the flap, this implies that 
there was some wash-out present in the model. This could be reduced in the future by 
increasing the torsional stiffness of the trailing-edge flap or attaching the actuation 
mechanism to a more central location on the flap instead of at the inboard end. Regardless of 
this effect, the measured actuation performance met and exceeded the goal of ±10 degrees 
dynamically. Nearly 10 degrees can be maintained for up to 30 Hz at only 14 psi PAM 
operating pressure, whereas nearly 18 degrees can be maintained for up to 35 Hz when 
driving the PAM actuators with 28 psi. Recall that this test case is a reduced load from the 
expected condition, so the PAM input pressures had to be reduced, as well. Based on all of 
these results, it can be stated that PAM actuation systems have clearly demonstrated their 
high performance capabilities for aerospace applications. 

5. Conclusion 

This research has developed and tested a series of innovative trailing-edge flap actuation 
systems that exploit antagonistic configurations of Pneumatic Artificial Muscles (PAMs) to 
generate bi-directional flap deflections. The systems were designed and built for 
experimental evaluation on the bench-top under simulated aerodynamic loadings with 
spring mechanisms and in the wind tunnel under actual aerodynamic conditions up to the 
maximum speed (Mach 0.3) of the Glenn L. Martin wind tunnel at the University of 
Maryland. Results showed that the flap deflection range produced was attractive to various 
flight control regimes, including flight control, vibration control, and even noise control. The 
key conclusion of this work is that PAM actuation systems have demonstrated the ability to 
dynamically control large flap deflections over a wide bandwidth in these varying control 
regimes and offer an attractive solution to aerodynamic control applications. 
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1. Introduction 



The tremendous destruction caused by recent typhoons in Japan has caused a substantial 
upsurge in interest in the subject of global warming among news media and the wider 
public. There are concerns that global climate change may have played a significant role in 
these events. Some believe that global warming is responsible for an increase in the 
frequency of destructive natural events. Typhoons cause the destruction of tiles on the 
rooftops of Japanese residences. The wind load on a roofing element is created by the 
difference between the external and internal pressures. The net wind load is, in general, 
determined by the building flow field, wind gustiness, and the element flow field (Peterka et 
al, 1997; Cermak, 1998). Although these parameters directly influence the external pressure 
distribution on a roofing element, the development of internal pressure, which indirectly 
depends on these parameters, is governed by a dynamic response that varies according to 
different roofing elements. The pressure distribution on an external roof surface and internal 
pressure distribution have been determined in numerous studies (Hazelwood, 1980; Ginger, 
2001). Element wind loading may differ significantly from the load derived from the 
external pressure distribution. Internal pressure is governed by the wind permeability of the 
surface, which is determined by openings, such as gaps between tiles or venting devices, 
and by the equilibrating resistance through and underneath a wind permeable surface 
(Kramer etal, 1979). 




Fig. 1. Japanese residence and roof tiles 
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Flow-induced vibration of roof tiles usually appears just before they are scattered. The flow- 
induced vibration (aeroelastic instability) of structures is an important phenomenon for the 
following two reasons: (1) strong lateral self-excited oscillations can develop at a certain 
wind velocity (onset velocity) as a result of the lateral aerodynamic force component and (2) 
these vibrations have a tendency to affect the behavior of the structure prior to the onset 
velocity because they produce negative aerodynamic damping that can considerably reduce 
the total damping available to the structure (Naudascher et al., 1993). However, the flow- 
induced vibration of roof tiles prior to scattering has been given very little attention. This 
study investigates the nature and source of the vibrating and scattering behavior of the roof 
tiles in order to provide better insight into this mechanism. This paper presents the first 
results of studies on the wind-inducing mechanism in roof tiles, which are widely used for 
roofing Japanese wooden dwellings (Fig. 1). 
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Fig. 2. Outline of the research 

Using wind tunnel tests, an experimental study was conducted to explain the behavior of 
roof tile vibration and the primary factors that affect their scattering. The results indicate 
that the vibration mechanism behaves in a manner that is consistent with that of a self- 
excited system, and the surface flow creates reasonable up-lifting moments only when the 
wind direction is roughly perpendicular to that of the eaves (Fig. 2). 

Nomenclature 

pitch angle (degree) 

(/) flow angle (degree) 

LI upstream flow velocity (m/s) 

X streamwise coordinate 

y transverse coordinate 

Z coordinate perpendicular to the surface of a roof tile 



2. Test facility and analysis procedure 

Fig. 3 illustrates the general layout of the apparatus used in this experiment. The 
experiments were conducted in an open-circuit wind tunnel that was driven by an axial 
flow fan. The nozzle of the wind tunnel had a 500 mm x 1,300 mm cross section. The 
maximum velocity of flow from the nozzle was approximately 50.0 m/s. The representative 
wind velocity was measured by a hot-wire anemometer and a linearizer on the exit nozzle of 
the wind tunnel. Approximately 10.0% of the flow's streamwise turbulence intensity was 
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generated by the grids. The spatial characteristics of air jet were checked for uniformity in 
wind speed and turbulence to ensure that all tiles were exposed to a near uniform air flow. 
The turbulence intensity of the flow condition is of the same order as the turbulence 
intensity experienced in practice. 
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Fig. 3. Experimental apparatus 

25 roof tiles were set up in 5 rows x 5 columns on a pitched roof in the downstream flow of 
a wind tunnel (Fig. 3). The roof tiles were tested by the air flow, which barely covered the 
entire exposed area of the tiles. They were made of clay, and each weighed approximately 
2.8 kg. The tiled pitched roof was fitted similar to a real roof arrangement with a plenum 
underneath the tiles, which acts as a roof cavity. This plenum was sealed with a clay pad. 
The internal pressure in this plenum was monitored and regulated by a pressure transducer 
placed underneath the tiles. The vibrations of the roof tiles were measured by a laser 
Doppler vibrometer (LDV, OMETRON VSIOOO) and an accelerometer (ONO SOKKI NP- 
3560, Fig. 4 (a)), and the normal natural frequencies of the roof tiles were analyzed using an 
impulse force hammer test. The vibration velocity could be measured up to 1,000 mm/s by a 
1 mW LDV, and the range of the vibrational frequency was from to 50 kHz. One roof tile 
was equipped with an accelerometer (Fig. 4 (b)). The accelerometer was used to measure the 
dynamic behavior of the tiles in three directions, X-, Y-, and Z under a no-flow condition. 
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(a) Impulse force hammer b) Frequency response function and coherence function 

Fig. 5. Frequency response function and coherence function of a roof tile generated by an 
impulse force hammer test 

and weighed approximately 5.0 g. The experimental measurement of the vibration 
frequencies for tiles was performed with the accelerometer. However, the vibration 
frequencies identified by the LDV were limited to small-amplitude modes. In this study, the 
accelerometer and LDV were used to determine the resonant frequencies of roof tiles that 
were and were not bolted to the roof bed. 

An impact hammer with a force transducer was used to excite the tiles under no-flow 
conditions (Fig. 5 (a)). Two signal conditioners were used to provide power to the 
accelerometer and the force transducer, whose spectral analyses were performed using a fast 
Fourier transform (FFT) spectrum analyzer (ONO SOKKI DS-2100 4CH). The sampling 
frequency was 5,120 Hz over a frequency range of - 2.5 kHz; 1,024 data points were sampled 
per spectrum. Unless otherwise stated, 64 spectra were averaged for each measurement. The 
frequency resolution of the spectra was 5 Hz. In order to analyze acceleration in a direction 
perpendicular to the surface of a roof tile, the time taken by the acceleration signal was 
recorded using the FFT analyzer. Two accelerometers were used simultaneously. Roof tiles 
that showed significant vibrations at any velocity, found from a series of experiments using 
accelerometers, were attached to two neighboring roof tiles on a model roof. 
The dynamic instability of the structure under excitation was also visually investigated. Large 
amplitude vibrations and the scattering of roof tiles were observed by a high-speed video 
camera (PHOTRON FASTCAM-PCI 2KC). The images were acquired at 2,000- frames per 
second, at a resolution of 512 pixels x 480 pixels per full frame. A hot-wire anemometer and a 
linearizer were used to measure the turbulence intensity of surface flow over the roof tiles. 



3. Results and discussion 

3.1 Impulse force hammer test for roof tiles 

Fig. 5 (b) shows the frequency response function curve and coherence function curve of roof 
tiles measured using an impact hammer with a force transducer. One of the resonant 
frequencies obtained by the accelerometer was 478 Hz. As stated in the next section, the 
measured frequencies obtained using the wind tunnel test are nearly consistent with the 
resonant frequencies obtained by the excitation analysis of the impulse force hammer test. 
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The value of the input excitation level is set to be approximately constant for the excitation 
analysis. However, the flow-induced excitation level is amplified and a higher level should 
be provided to obtain vibration measurements. On the other hand, the variation in the 
measured values of resonant frequencies for the accelerometer measurement and excitation 
analysis may be attributed to the added weight of the accelerometer in this experimental 
technique. In order to eliminate the effect of the added weight of the accelerometer on the 
resonant frequencies of the roof tile, the corresponding frequency response curve of this roof 
tile was obtained using the LDV. The peak values of this frequency response curve were 
compared with those obtained using the accelerometer method. It was found that the results 
of resonant frequencies measured using LDV and those using the accelerometer agreed 
satisfactorily. 

3.2 Acceleration measurements of roof tile 

In the measurement and analysis of roof tile vibration and its acceleration, the pitch of the 
roof 6 was set at 19 degrees, 24 degrees, and 29 degrees and the flow angle ^ was set at 
degrees. The wind velocity was gradually increased from to 50.0 m/s or until scattering of 
the tiles occurred. The signals from the accelerometer s were recorded to be analyzed later 
using a personal computer. 

The slope angle of the roof was changed, and the effects of fluttering in the last stage of roof 
tile scattering were examined (Figs. 6-8). The small-amplitude vibration of the model roof 
tiles appeared in a low-velocity flow at the maximum pitch angle of 29 degrees, while the 
model roof tiles showed fluttering when the wind velocity reached approximately 38 m/s. 
They were more buffeted at the pitch angle of 24 degrees than at the pitch angle of 29 
degrees, and then fluttered when the wind velocity reached approximately 40 m/s. The 
model roof tiles did not flutter at the minimum pitch angle of 19 degrees, and they were 
buffeted at a higher wind velocity than that at other pitch angles. They did not flutter at 
pitch angles of 24 and 29 degrees because of the weight of the neighboring roof tiles and 
bolts. The fluttering of the model roof tiles appeared with relatively large-amplitude 
vibrations, and it was regarded as fluttering when the roof tile was completely lifted from 
the roofing board and the board was exposed. 
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Fig. 6. Effect of slope angle of roof on vibration of roof tiles at ^ =29 degrees, LI = 39.0 m/s 
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Fig. 7. Effect of slope angle of roof on vibration of roof tiles at i9 =24 degrees, LI = 38.5 m/s 
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Fig. 8. Effect of slope angle of roof on vibration of roof tiles at i9 =19 degrees, LI = 39.9 m/s 




Fig. 9. Observation of the flow on the surface of the roof tile by the oil film method 
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The occurrence of fluttering was considered as one of the last stages of roof tile scattering. 
Fluttering did not occur at the lower pitch angle and, as a result, the model roof tiles resisted 
vibration and scattering at even higher velocities, whereas fluttering occurred at a higher pitch 
angle because the model roof tiles were often buffeted and scattered by lower critical 
velocities. In addition, the oil-film method was used to observe the flow pattern on the surface 
of the every tile in the model roof (Fig. 9). Separation regions appeared over the surface of the 
roof tiles as the flow angle was gradually increased. The wind flow was along the surface of 
the roof at a lower pitch angle, whereas the flow was split directly by the surface and the edge 
of the roof tile at a higher pitch angle. The wind flow was toward the edge of the roof, and it 
was found that the flow became very turbulent locally over the edge and formed the 
separation regions. The appearance of the separation regions was considered to indicate a 
significant fall in the external pressure and a rise in the internal pressure, causing the fluttering 
of the roof tiles. The clattering of the model roof tiles was recognized in the Y-axis direction at 
a pitch angle of 29 degrees. It was found that this may be caused by the decrease of the critical 
velocity, which lifted up the tiles and caused them to flutter. 

Figs. 6-8 show the results for several pitch angles. The largest number of roof tiles were 
buffeted and fluttered at the maximum pitch angle of 29 degrees (Fig. 6 (a)), resulting in 
three roof tiles being scattered. The roof tiles were not scattered at a pitch angle of 24 
degrees (Fig. 7 (a)), whereas they were both buffeted and scattered at a pitch angle of 29 
degrees. In this case, the roof tiles attached with the accelerometers and the eight 
neighboring roof tiles fluttered. At the minimum pitch angle of 19 degrees (Fig. 8 (a)), the 
roof tiles did not flutter even at the maximum velocity of 40 m/s, but a few roof tiles were 
buffeted. 

Similarly, the signals from the accelerometers revealed fluttering at pitch angles of 24 and 
29 degrees. The results obtained by the FFT analysis of the accelerometer signals at pitch 
angles of 24 and 29 degrees are shown in Figs. 7 (b) and 6 (b), respectively. The small- 
amplitude vibrations were recognized using accelerometers 1 and 2 at a low velocity at 
these pitch angles. However, a 50 m/s^ vibration was recognized momentarily when the 
wind velocity was increased. The small-amplitude vibrations appeared at a low velocity, but 
gradually increased at a higher velocity (Peter ka et al., 1997; Cermak, 1998). Moreover, 
acceleration and amplitude of the vibrations are related to the critical velocity at which 
fluttering occurs for a constant pitch angle of the roof. At a pitch angle of 24 degrees, the 
values of acceleration and amplitude increased when the wind velocity reached 40 m/s (Fig. 
7 (b)). The acceleration and amplitude observed at a relatively high velocity for a pitch angle 
of 24 degrees also appeared at a low velocity for a pitch angle of 29 degrees (Fig. 6 (b)). 
The maximum acceleration value prior to fluttering was found at pitch angles of 24 and 29 
degrees. Moreover, it was found that the acceleration decreased when the roof tiles fluttered 
at both pitch angles. This was found to be caused by the balancing of internal pressure in the 
space between the attic side and the roofing board by the external pressure of the flow over 
the roof tiles. On the other hand, it was occasionally observed that the neighboring roof tiles 
touched each other and unexpected acceleration signals were found because of contact 
between neighboring roof tiles during vibration. The effects of the roofs pitch angle on their 
scattering were recognized by analyzing the acceleration signals. 

In a series of wind tunnel tests, the pitch angles of the roof were changed and two 
accelerometers were attached to neighboring roof tiles in order to detect and analyze 
acceleration signals. Fig. 10 shows an example of the acceleration signals that were 
frequently found whenever small-amplitude vibrations occurred. It shows the acceleration 
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signals and the behavior of the 25 model roof tiles at a pitch angle of 24 degrees and wind 
flow of 40 m/s. In this experiment, it was also recognized that the tiles locally arranged at 
the back and right side of the roof were often buffeted and scattered by strong winds. In a 
typical construction method, tiles are piled up and laid on the upper and lower side of a roof 
by their own weight. Therefore, the roof tiles locally arranged at the back and right side of 
the roof are held by the relatively lighter weight of the neighboring roof tiles. The wave- 
forms shown for case A in Fig. 10 (a) indicate that the signals are out of synchronization by a 
half cycle from both accelerometers. The wave-forms shown for case B in Fig. 10 (b) indicate 
the nearly synchronized signals from both accelerometers. Fig.lO shows the time history of 
the acceleration signals caused by the wind effect. These results show that it was found that 
the accelerations and time histories of the two roof tiles differ from each other. The forces 
affected both roof tiles simultaneously. On the other hand, even when a force did not act 
directly on a roof tile, it affected the other roof tile. In this case, the roof tile equipped with 
accelerometer ® was held by the roof tile equipped with accelerometer (D and the force of 
the roof tile equipped with accelerometer (D was added, although the force of the roof tile 
equipped with accelerometer ® did not act directly on it. The force of the tile equipped with 
accelerometer ® was added to the roof tile equipped with accelerometer (2) although the 
force of the roof tile equipped with accelerometer (2) did not act directly on it (case A). It was 
observed during a series of experiments that the roof tiles were buffeted almost 
simultaneously when vibration occurred (case B). 
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Fig. 10. Acceleration signals from accelerometers 1 and 2 at ^ =24 degrees, LI = 40.0 m/s 
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These results suggest that the roof tile held by the other roof tile was buffeted and then 
produced the reaction force on the neighboring roof tile. This behavior may be same as the 
"wave motion of roof tiles'' phenomenon described by Ginger (2001), which is often 
reported to cause construction damage. 

Hazelwood (1980b) described the lifting mechanism of a tile by a moment turning the tile 
upward around the pivoting point on the batten. The moment consists of a lifting force and 
two force couples caused by the external and internal pressure distributions, respectively. 
Fig. 11, which shows snapshots of a high-speed video camera picture, demonstrates this 
lifting mechanism for wind direction perpendicular to the eave. For a local flow direction 
perpendicular to the ridge, the internal pressure in the space between the tiles and underlay 
may become positive because of stagnation. As a results, the net wind load increases 
because of the sealing effect of the underlay. However, pressure equilibration in the gable 
roof is prevented, leading to a much lower net wind load for the leeward roof tiles. 
Aerodynamically favorable tiles should have a shape that prevents stagnation at the 
overlaps. The permeability at the overlapping gaps parallel to the ridge should be small and 
the permeability at the interlocking gaps perpendicular to the ridge, where suction occurs 
due to the element flow field, should be high (Hazelwood, 1980a). 




Fig. 11. Lifting up of roof tiles due to wind action as recorded by a high-speed video camera 
atLZ = 40.0m/s 



3.3 Transverse vibration of roof tiles 

Fig. 12 shows an example of the typical turbulence spectrum obtained by the hot-wire 
anemometer for 0=0 degrees, ^ = 90 degrees, and LI = 40 m/s. With a turbulence level 
of surface flow close to the roof tiles, the tiles exhibited only the typical turbulence- 
buffeting response within the intermediate ranges of the angle of incidence. The Reynolds 
number during the experiment was so high that the edge separation was turbulent. The 
sources of vibration are the front and side edge vortices (Fig. 13). The vibration 
amplitudes increased progressively with increasing velocity, which indicates a typical 
buffeting response. 
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(c) Vibrational acceleration and turbulence power spectrum. 

Fig. 12. Vibrational acceleration and turbulence power spectrum of roof tiles for 9=0 
degrees, ^ = 90 degrees, and LI = 40.0 m/s 
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Fig. 13. Transverse vibration of roof tile 

The local flow due to the outer shape of a surface element is of importance if the element is 
located in an area with attached flow, such as on the windward surface of a pitched roof. 
The gaps between the tiles may be exposed to local stagnation and/ or suction depending on 
the shape of the tiles. If suction prevails, the internal pressure is decreased and the opposite 
takes place for predominating stagnation. For 6> = 30 degrees, a front edge vortex with its 
axis parallel to the ridge is formed, causing significantly higher negative pressure 
coefficients (Ginger, 2001). It is observed by the surface oil-flow visualization method that 
reattachment takes place upstream of the ridge and the flow is completely separated at the 
leeward roof area. If the roof pitch is increased, the vortex on the windward side decreases 
in size and reattachment takes place much closer to the eave. In the region of flow 
bifurcation, the pressure coefficient becomes positive (Peterka et al., 1997). 
However, if the external pressure distribution is changed because of the shape of the 
element, the internal pressure can be affected significantly. In particular, for the local flow 
direction perpendicular to the ridge of a tiled roof, the flow is stagnated at the overlaps of 
the tiles. The stagnation pressure increases because of the step formed by overlapping tiles 
and leads to an increase in the internal pressure if the permeability of the overlap gaps is 
sufficient. This value depends on the shape of the front and side edges of the tile, i.e., square 
or round, and the level of free-stream turbulence; the larger the value of free-stream 
turbulence, the larger is the critical value of incidence. Because the pressure distribution on 
the roof is strongly influenced by the turbulence of the oncoming flow, this turbulence will 
also affect the net loading on roof elements. 

If a roof tile is inclined with respect to the free stream, the flow will separate from one side 
as soon as the angle of incidence exceeds a critical value. Visualization using the surface oil 
flow method shows that the vortex cones caused by the yawing flow separation at the 
leading edges result in the highest negative pressure coefficients close to the windward 
gable and the windward eaves. If the roof pitch is increased, the vortex cones decrease in 
strength. In regions of separated flow, the external pressure distribution on a tiled surface 
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coincides with the pressure distribution on the roof surface, as described by Peterka et al. 
(1997). In regions of attached flow, however, the pressure distribution on a tile is influenced 
by the flow around the tile. A typical example for the change in the external pressure 
distribution due to the element flow field is shown in Hazelwood (1980a). 
The pressure distribution, indicating an acceleration region at the eave-facing end of the tile 
and a stagnation zone in front of the overlap of the tile in the upper row, results in an 
upward-lifting moment. The predominant geometric parameter for the pressure distribution 
is the tile thickness related to the non-overlapping length (Peterka et al., 1997). The 
fluctuations of the surface flow velocity caused by the instabilities of the flow field over the 
roof will change the pressure distribution and make the tiles clatter. When the wind load 
exceeds a certain value, the tiles are lifted up and the permeability of the roof surface 
increases rapidly. If this happens in a region with low external pressure, the wind load on 
the tiles will decrease. However, if lifting-up occurs because of surface flow action on the 
windward side, the stagnation effect will lead to an increase in the internal pressure and the 
up-lifting tile load. The internal pressure underneath the tiles affects the overall stability of 
the tiles and acts as the up-lifting tile load. 

The small-amplitude vibrations of the roof tiles appeared first, the amplitude grew 
gradually larger as the wind velocity increased, and then fluttering with large-amplitude 
vibrations occurred, finally followed by scattering. The vibrational frequency was identified 
by image analysis of the high-speed video camera to measure relatively high-amplitude 
vibration in fluttering, which is considered to be the direct cause of tile scattering. The roof 
tiles do not always oscillate with a fixed vibrational frequency. Because vibrations with 
several frequencies affected the tiles and showed complex behaviors, some oscillation 
patterns were chosen at random from the data to be analyzed further. It was found that the 
amplitudes of tile vibration were larger than that of their natural frequency, and the 
vibration frequencies were low (in the range of 10 - 20 Hz). 

The results obtained by the FFT analysis of the acceleration signals in the experiment in 
which fluttering occurred are shown in Fig. 14. The results show the oscillation of fluttering 
at a pitch angle of 24 degrees and a wind velocity of 40 m/s. The wind velocity was 
gradually increased from the start of the wind tunnel test to its maximum velocity, and the 
acceleration measurement and the video camera recording were then started 
simultaneously. The sampling time of the FFT analyzer was set at 2,048 points, the 
frequency resolution was set at 800 lines, and the frequency range was 0-5 kHz. Moreover, 
the peak frequency of approximately 470 Hz, which appeared just before tile scattering, was 
the natural frequency and was also recognized by FFT analysis. To minimize the effects of 
sampling time on the results of the FFT frequency analysis, the FFT frequency was analyzed 
using a sufficient sampling time. As a result the relatively high frequency, i.e., the natural 
frequency, as well as the relatively low frequencies were recognized. 

It was observed in the wind tunnel test that the bolted roof tiles were lifted up, damaged, 
and then scattered by the wind, and they induced further fluttering and clattering by lifting 
up their neighboring roof tiles. In other words, it is believed that the amplitude was the 
largest in one cycle of tile vibration and the largest energy was obtained at those moments. 
The force acting on the roof tile can be estimated by Newton's second law of motion. In the 
case of the measured acceleration of 11 m/s^ and the given mass of 2.8 kg, the force acting 
on the roof tile was 30.8 N. 
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Fig. 14. Vibrational acceleration power spectrum of roof tiles at ^ =24 degrees, 
LI = 40.0m/s 

The natural frequency of the roof tile was measured by the impulse force hammer test. The 
center of a roof tile hung from the ceiling was hit by the impulse hammer. The natural 
frequency of the tile was analyzed in terms of a frequency-response function and a 
coherence function. By analyzing the frequency-response function, the peak frequency was 
found to be 478 Hz. The coherence function was strongly correlated with the frequency- 
response function (Fig. 5 (b)). It was recognized that the dominant frequency, which 
occurred just before the scattering shown in Fig. 14, almost coincided with the natural 
frequency of the tiles that was found by the impulse force hammer test. The natural 
frequencies of the roof tile hung from the ceiling were found to be between 430 and 460 Hz. 
The peak frequency of the roof tile appeared just before scattering, as shown in Fig. 14. The 
roof tiles were arranged on the model roof in order to measure their vibrational frequency 
caused by the wind at the center of the opposite side of the roof. It was found that the 
measured frequency was different from the frequency of fluttering and the natural 
frequency of the tiles (Naudascher et al., 1993; Hazelwood, 1980b). 

These test results showed that the vibrational frequency of about 14 Hz almost coincided 
with the vibrational frequency that was obtained by analyzing the images of the high-speed 
video camera. On the other hand, the information of the acceleration and the results of the 
image were analyzed to specify the vibration occurring during fluttering. Low-frequency 
vibrations (10 - 20 Hz) were detected first (Fig. 14). Next, the significant peak amplitude of 
the natural frequency, which appeared just before fluttering, was also recognized. In other 
words, it is believed that the vibration at the relatively low frequency has a dominant effect 
on fluttering, and this natural frequency appears prior to fluttering because of the significant 
vibration at the relatively high natural frequency just before fluttering. Finally, the 
occurrence of vibration at the low frequency with a relatively large amplitude has the 
greatest effect on fluttering, and this mechanism can result in the lifting of the roof tiles. 
Hence, the dynamics of the roof tiles were due to the balance of their own weight, to which 
the external pressure was added by the fluid over the surface of the roof, and the internal 
pressure (i.e., the space between the roof tile and the roofing board). Because the external 
pressure and the internal pressure were changed, an unbalance of both pressures occurred, 
the tiles became unstable, and then fluttering occurred. It is believed that the relatively low- 
frequency vibrations have the greatest effect on scattering and can be the main factor that 
controls the behavior of the roof tiles. 
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4. Future research 

Strong winds not only result in tile scattering leading to damage of tiles, but also result in 
water leak damage. Experiments pertaining to water leaks can be broadly classified into 
pressure box- type experiments and blower/ water dispersion- type experiments. Pressure 
box-type experiments allow for the recreation of model wind pressures using devices for 
either increasing or decreasing pressure. Conversely, blower/ water dispersion- type 
experiments make use of devices consisting of blowers and water dispersion equipment, 
which allow experiments to be conducted in conditions very similar to the actual flow of 
wind and rain during stormy weather. However, these types have both advantages and 
disadvantages, and neither of these types is able to reproduce the actual conditions of both 
rain damage and the damage caused by heavy winds simultaneously. 

In future work, the authors will focuse their attention on vibrations which cause the 
preliminary phenomena eventually leading to the scattering of tiles due to the effects of the 
wind, and will seek to understand the mechanism of these vibrations. Consequently, the 
authors will be able to connect together the mechanism and the preliminary phenomena of 
the occurrence of tile vibration induced by fluid flow. In accordance with these results, in 
future work the authors will go beyond the conventional understanding of water leakage 
amounts, aiming to establish appropriate experimental methods and to clarify the 
mechanism underlying the occurrence of water leak phenomena. The authors intend to 
investigate the previously unknown influence exerted by tile vibrations on water leaks. The 
ultimate goals are to provide a conclusive understanding of the effects of wind and to 
provide suggestions for possible improvement and redesign of roof tiles (Fig. 15). 

5. Conclusions 

An experimental study was conducted using wind tunnel tests in order to explain the 
behavior of roof tile vibration and the primary factors that affect scattering. The results are 
summarized as follows. 

1. The basic mechanism that can lead to flow-induced vibrations of roof tiles is similar to 
that of the so-called fluttering instability, which appears as self-excited oscillations in 
the natural mode of a structure at a certain critical flow speed. The oscillating 
frequencies are related to the natural frequencies of vibration. 

2. Surface flow is only important on the windward side of a roof and creates reasonable 
up-lifting moments only for wind directions roughly perpendicular to the eaves. 

3. The effects of a roofs pitch angle on the fluttering of roof tiles were confirmed 
by analyzing acceleration information as the pitch angle was increased; the absolute 
value of acceleration and the amplitude also increased with increasing pitch angle. 

4. The ''wave motion of roof tiles'' appeared just before scattering was observed, and the 
forces acting on two neighboring roof tiles were found to be either synchronized or out 
of phase. 

5. Low-frequency vibrations, which have the greatest effect on scattering, were identified 
by a high-speed video camera, and the major factor that controls the behavior of the 
roof tiles was found to be the balance between the external pressure and the internal 
pressure. 
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